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A  subsequent  series  of  test  runs  performed  on  the  SUBSONIC  program, 
indicated  that  the  results  were  not  reliable  over  a  large  range  of  flight  conditions.. 
Primarily,  this  was  due  to  a  lack  of  refinement  in  the  original  analysis  pertaining  to 
the  inconsistency  incolving  y  along  the  reference  line  (page  20).  In  addition,  an 
instability  condition  affecting  the  system  of  pertubation  equations  when  a  large  number 
of  iteration  points  were  prescribed,  was  not  detected  initially. 

In  the  course  of  solving  these  and  other  problems,  an  exhaustive  series  of 
tests  were  conducted.  Both  the  analysis  and  the  coding  were  greatly  refined;  indeed, 
a  reiteration  procedure  was  added  which  greatly  enhances  the  accuracy  of  the  final 
solution,  at  a  very  small  additional  cost  in  machine  time.  A  scaling  technique  was 
introduced  which  not  only  improved  the  solution  in  the  elliptic  region,  but  insured 
that  a  subsequent  execution  of  the  supersonic  flow  field  would  have  accurate  initial 
conditions,  consistent  with  the  subsonic  flow  properties.  The  precision  of  the 
calculated  deviations  along  the  body  profile  was  considerable  improved  by  refining 
the  integration  along  the  reference  line. 

These,  and  many  other  additional  features,  have  yielded  a  computer 
program  which  will  generate  accurate  results  with  a  high  degree  of  reliability  at 
relatively  low  cost. 
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Page  3: 

Last  sentence,  second  paragraph,  replace  with 

"  The  coefficient,  a^  is  then  modified  parametrically  to  overcome 
this  situation." 

Last  sentence,  third  paragraph,  replace  with 

"  The  reference  line  is  located;  then  the  "basic  polynominal"  is  scaled 
so  that  the  body  point  on  the  reference  line  falls  on  the  prescribed  profile.  See  page 
19  for  details. " 

Page  4: 

Last  sentence,  third  paragraph,  add 

This  "final"  polynominal  (and  reference  line)  is  then  scaled. 

6.  To  refine  this  solution,  a  second  set  of  sweeps  is  executed  (steps  2-5) 
using  the  first  "final"  polynominal  as  the  second  "basic  polynominal"  and  utilizing 

a  new  set  of  "iteration  points". 

The  .last  paragraph  starts  as.  follows: 

7.  The  Tree  stream  conditions,.... 

Page  18: 


Delete  last  4  lines 


ii 


Page  20: 

Replace  the  last  2  sentences  with, 

"These  inconsistencies  are  minimized  by  assigning  as  the  y  used  in  the 
Subsonic  region,  the  average  of  the  values  at  the  body  point  and  behind  the  detached 
shock  at'  the  reference  line.  In  addition,  the  coordinates  of  the  body  point  are  defined 
by  the  condition,  1”  =  0. 

Page  22: 


Page  24: 


Page  26 : 


Page  27 : 


Lines  9  and  1 6 ,  replace  with  Y^^, 


Delete  the  next-to-last  sentence  in  the  first  paragraph. 


The  exponents  on  lines  4  and  6  should  be  -  V2  (  not  +  V2  ). 


Replace  lines  6,7,  and  8,  as  follows: 


c.  =  C. 
m  in 


d.  =  D. 
m  in 


( T  -T  ) 
n  ns 


(y 


|i+  ij  ^1- 1) '  j 

^n'TJ 


l'KH’ 


V 


e.  =  -  2  ( A.  U  +  R  V  )  -  (  A.  B.  -  E.  )  ( T  -  T  J 

m  m  ns  m  ns  m  clTns  m  '  n  ns 


Ill 


Page  28: 

Delete  the  first  4  lines  and  replace  equ.  (4)  with 


3_U 

a  y 


U(n+1)  ~  U  (  n- 1  ) 


S  1  £  n  :£  p-1 


Page  30: 


Add  to  right  of  line  5, 

;  a,  j8  given  on  page  26. 
Add  at  bottom  of  page, 


Then, 


(a  u| 

’  Un+1  -  Un-1  " 

\a  y  J 

n 

.  yn+l  '  yn-l  J 

la  v| 

[  Vn+1-Vn-1  ] 

y 

n 

.  yn+l  '  yn-l  . 

Page  31 : 

g 

Delete  the  term,  +  ag  y  ,  in  line  2 

Replace  the  symbol,  6  »  with  6  >  in  line  9  and  in  sketch. 


Page  46  : 


Line  6,  add  the  term, 


e'/Tcp 


e 


(4) 


IV 


Page  68: 

Line  12,  insert,  ....  pages  37  through  41. 

Page  72: 

Equ.  (16)  :  x  (  2  0  )  =  |  (  -  xq  )  +  xQ 


Page  75: 

Add  this  sentence  at  the  bottom: 

It  is  necessary  to  "right- adjust"  data  words  within  the  indicated  field; 
that  is,  the  word  must  be  shifted  to  the  extreme  right  of  the  field. 

Page  79: 

Delete  the  last  sentence  of  paragraph  4. 

Page  85: 

Replace  the  2  with  a  3  on  next  to  last  line. 

Page  86: 


Add  this  sentence  at  the  bottom: 

Note  that  the  Subsonic  Program,  coupled  with  this  option  of  the  Supersonic 
Program,  is  exactly  equivalent  to  the  Flow  Field  program. 

Page  98: 


Lines  4,  5,  and  6  are  to  be  replaced: 

P  Dimensionless  pressure  = 


_ E _ , 

p  w  1 

oo  oo 


R  =  P  =  6  -f- 
^00 


Pressure,  lbs/sq.  ft. 

Dimensionless  density, 

where  y  =  */2  [y  + 
'  L  shock 


6  = 


r  -  1 

y  +  1 


line 


2a 


The  range  of  applicability  of  these  programs  is  limited  only  by  the 

numerical  fits  supplied  by  ABMA,  which  express  the  entropy  and  density  in  terms 

of  the  known  enthalpy  and  pressure.  Those  flight  conditions  of  velocity  and  altitude 

o 

at  which  the  stagnation  temperature  is  less  than  7000  K  and  stagnation  pressure  less 
than  10  atmospheres,  lie  within  the  scope  of  these  programs. 
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Note  that  the  preceding  analysis  is  based  upon  the  assumptions  of 
linear  mass  flow  variation  between  shock  and  body  at  the  axis,  a  parabolic 
variation  across  the  shock  layer,  r  s  ,  and  a  cubic  variation  at  the  down¬ 
stream  section,  c  d.  These  assumptions  are  consistent  with  the  physics; 
as  one  proceeds  around  the  "shoulder",  in  the  downstream  direction  this 
variation  assumes  still  higher  orders. 

After  sweeping  through  the  transonic  region  using  this  "basic"  poly¬ 
nomial,  the  properties  at  the  point  on  the  body  at  the  reference  line,  are 
compared  with  the  prescribed  body  profile.  In  general,  the  calculated  values 
of  x,  y,  and  6  at  this  point  (  which  we  will  call  x,  y,  0  )  will  not  be  consistent 

with  the  prescribed  profile,  We  can  compute  the  coordinates  of  the  point  on 

—  Yr 

the  prescribed  body  (  xg,  yg)  ,  where  Qg=  9,  Then  define  K  =  — — — .  Now 

y 

this  "basic"  shock  polynomial  may  be  modified  so  that  this  calculated  point 

on  the  body  agrees  with  the  prescribed  body  profile.  If  we  now  denote  the 

original  coefficients  by  an,  n=  0,2,4  instead  of  an,  we  can  define  the 

modified  polynomial  as  x=  F(  y)  =  )  anyn  and  relate  them  as  follows; 

n=  Oj  2,  4 


x  =  F(  y)  = 


Ka 


o+  (^b'Kx: 


»1  I 


( 1-n)  _  n 

+  i_  K  anv 

n=  2,4 


a0=Kao  +(xBrKx):  an  = 


K 


(  1-n) 


^n  ’  n=  ^ 


Thus, 
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SUBSONIC  ANALYSIS 

At  the  completion  of  tke  transonic  region,  the  reference  line  is  defined, 
and  the  program  then  transforms  the  properties  of  all  points  on  the  reference 
line,  from  the  physical  x-y  plane  to  the  T -y  plane.  The  properties  along  the 
reference  line  required  by  the  Subsonic  analysis  are  U,  V,  T  and  y. 


U  = 


W 


W 


cos  0 


oo 


00 


sin  0 


y  =  y 


The  values  of 


W 


at  each  point  are  calculated  by  imposing  the  condition 


that  the  values  of  pressure  along  the  reference  line  calculated  by  the  method 

of  characteristics  in  the  transonic  region,  and  by  the  subsonic  analysis  in  the 

elliptic  region,  be  consistent.  The  dimensionless  pressure,  P,  used  in  the 

subsonic  analysis  is  defined  by  the  expression,  P  =  .  P  ;  calculate 

Poo  WooT 

P  where  p  is  the  value  of  pressure  in  psf  at  points  along  the  reference  line, 
calculated  by  the  method  of  characteristics.  By  transposing  the  expression 
at  the  bottom  of  page  25 

,1/2 


Vu 


2+v2 


W 


c  -  (  1  +  6)  f  f 


2-6 

1+6 


The  above  approach  is  necessary  since  the  value  of  y  varies  from  point 
to  point  along  the  reference  line  as  computed  by  the  method  of  characteristics, 
while  y  (  and  6)  is  held  constant  in  the  subsonic  analysis. 
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T  is  found  by  integrating  along  the  reference  line  from  the  shock  (T«  1) 
towards  the  body  according  to  the  following  formula.  Refer  to  Fig.  II . 


t  -  nr  & 

1 ‘ns 

Yn+Yns  } 

— 2 '  yn  yns^ 

+ 

(R‘U)n  +  (R-U)ns_ 

(  Yn_Yns)  "(xn"xne) 

(R-V)„+(R.V)IIS 

r  . 

where 


R„ 


Y  =  1/2  (Yb  +  Y0) 


Yb  is  the  value  at  the  body  on  the  reference  line 

Y0  is  the  value  behind  shock  at  intersection  with  reference  line 


For  the  perturbation  technique  to  be  valid,  the  integration  of  T  along  the 
reference  line  must  be  extremely  accurate  and  it  is  necessary  that  all  parameters 


used  in  the  integration  must  be  self-consistent. 


Thus  R  is  set  equal  to 


The  integration  is  performed  using  the  subsonic  mesh  points  as  pivotal  points  and 
iterating  on  Tn  at  each  step  to  ensure  maximum  precision.  Since  the  integrated 
value  of  T  at;  the  body  will  not,  in  general,  be  zero,  it  is  necessary  to  adjust  the 
position  (and  properties)  of  the  terminal  (body)  point  of  the  reference  line.  Thi6 
last  step  is  iterated  until  the  coordinates  of  the  body  point  on  the  reference  line  is 
consistent  within  very  close  tolerance,  with  an  integrated  value  of  T  =  0. 
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If  a  total  of  J  points  are  prescribed  as  iteration  points  (  excluding  the 
stagnation  point  -  J  =  3  in  above  sketch  ),  the  program  executes  J  additional  sweeps 
through  the  Transonic-Subsonic  regions.  For  each  sweep,  one  coefficient  of  the 
"basic"  shock  polynominal  is  perturbed  and  the  resulting  deviations  calculated.  Thus 
for  the  jth  sweep,  j=l,  2,  . , .  . ,  J,  the  detached  shock. polynominal  is  of  the  form 
F(y)  =  aQ  +  a2  y2  +  a^  y4  +  Aa2j  y  ^  ,  and  the  deviations  are  prescribed  as 

*U)  •  _  ,  , 

0.  f  i  — 
l 

After  J  sweeps  have  been  completed,  a  system  of  simultcineous  equations 
may  be  written  in  the  unknowns,  C  ,  as  follows:  =  6^-  {  See  sketch  ) 

cx  (  6  j1-  $  +  C,  (  ip-  ip)  + . +  Cj  (  fif  \  =  -  6i(o)  N 

q  (  62  4°)  +  C2  (  6 (p-  6{p)  + . +  Cj  (  ip  -  6 (p)  =  -  p  o) 


V  C1  (fi 

J 


(l  -  6(°j  +  c2  ( ip-  i0)' 


)  + 


or 


I 

m=l 


r  i  *(°)  \  _  *  (°)  ,9  t 

5  -0  )  =  -  o  ,  n  =  1,2, ...,J 

m  n  n  n 


where 


and 


The  final  polynominal  is  then  written 

F  (  y  )  =  aQ  +  S2  y2  +  y4  +  ....  +  a2j  y 


2J 


&2-  =  a^.  +  C.  .  Aa2.  ;  A  a-2-  is  the  jth  perturbation 

J  J  J  J  J 

j  =  1,2,3, . ,J 

J 


=  a 

o  o 


<o> 


+  L 

j=i 


5°’  -6  (0) 
o  o 


A  final  sweep  through  the  Transonic-Subsonic  region  uses  F  (y)  to  describe 
the  detached  shock.  The  points  along  the  reference  line  are  scaled  (see  below)  and 
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written  on  binary  tape  B  3  along  with  the  free  stream  conditions  and  the  body  equations, 

and  the  properties  within  the  Subsonic  region  are  printed  as  output. 

The  value,  K  (see  page  19),  is  calculated  for  this  final  reference  line, 

and  the  position  of  each  point  is  adjusted  as  follows: 

x  =Kx;1,  +  (x„-Kx);y  =  K  y  .  , 

new  old  B  'new  •'old 

The  "final"  shock  polynominal  is  modified  as  indicated  on  page  19,  and 


the  final  pass  through  the  elliptic  region  is  then  completed, 
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1 

CARD 

COL, 

DATA 

FORMAT 

L 

*- 

1 

4,  5, ....  A- 1 

.  2-10 

a 

F 

| 

12-20 

b 

F 

22-30 

c 

F 

32-40 

d 

F 

L 

42-50 

e 

F 

r 

52-60 

Xt 

F 

62-70 

a 

F 

A 

1-72 

Blank 

A  +  1 

1 

Stagnation  Point 

Code  =  1 

I 

i 

2-10 

xo 

F 

1. 

11 

Point  1  Code  = 

1  or  2 

I 

T 

1 

12-20 

X1 

F 

rr 

21 

Point  2  Code  = 

1  or  2 

I 

- 

22-30 

x2 

F 

n 

1 

31 

Point  3  Code  = 

1  or  2 

I 

* 

32-40 

X3 

F 

V.- 

41 

Point  4  Code  = 

1  or  2 

I 

1 

42-50 

xA 

F 

A  - 

4 

51 

Point  5  Code  = 

1  or  2 

I 

52-60 

X5 

F 

.  _ 

61 

Point  6  Code  = 

1  or  2 

I 

62-70 

X6 

F 

- 

A+2,  A+3 

Identical  in  format  to  card  A  + 

1 

(if  required) 


Pertains  to  points  7  through  20,  in  sequence 
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CARD 

COL. 

DATA  FORMAT 

B 

2-10 

y0 

F 

12-20 

Y1 

F 

22-30 

*2  ‘  ' 

F 

32-40 

y3 

F 

42-50 

y4 

F 

52-60 

y5 

F 

62-70 

y6 

F 

B+l.B+2 

Identical  in  format  to  card  B. 

(if  required) 

Pertains  to  points  7  through  20,  in  sequence 

C 

1,2 

Number  of  intervals  of  mesh  spacing,  A 

I 

3-10 

ATi 

F 

11,  12 

Number  of  intervals  of  mesh  spacing,  AT., 

I 

13-20 

AT2 

F 

21,  22 

Number  of  intervals  of  mesh  spacing, 

I 

23-30 

AT3 

F 

41, 42 

Number  of  intervals  of  mesh  spacing,  Ay^ 

I 

43-50 

Ayi 

F 

51-52 

Number  of  intervals  of  mesh  spacing,  A 

I 

53-60 

Ay2 

F 

61-62 

Number  of  intervals  of  mesh  spacing,  A 

I 

63-70 

Ay3 

F 

C  +  1 

1-10 

The  x-coordinates  of  the  "'.iteration  points” 

F 

11-20 

along  the  prescribed  profile  in  ascending  order 

• 

for  the  second  series  of  sweeps.  A  maximum 

• 

61-70 

of  7  points. 
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i 

8 

1 

I 

I 


CARD 

C  +  2 


COL. 

1-10 

11-20 

61-70 


DATA 


The  y-coordinates  corresponding  to 
the  previous  card 


FORMAT 

F 


,  80 
only.  The  corresponding  y-coordinates  are  on  cards  (s)  B  (  and  B+l  and  B+2  ). 

Note  that  the  stagnation  point  must  jjlways  be  prescribed  as  an  "iteration  point". 

To  assure  an  accurate  solution,  it  is  advisable  to  prescribe  at  least 
two  points  adjoining  the  assumed  sonic  point  on  the  body,  as  geometry  points.  For 
a  spherical  nose,  these  points  should  be  located  where  the  slope  of  the  body  is 
approximately  40°  and  50°  respectively.  Should  the  body  exhibit  a  more  rapid 
change  of  curvature  than  does  a  sphere,  in  the  neighborhood  of  the  sonic  point,  then 
additional  geometry  points  should  be  clustered  in  this  region. 

The  geometry  of  the  body  should  be  "normalized"  so  that  the  assumed 
sonic  point  has  the  coordinate,  y'J'  =  hg  £  1. 0  (see  sketch,  page  8). 

The  input  card,  C,  prescribes  the  T  -  y  mesh  used  in  the  subsonic  region. 
As  indicated,  up  to  3  different  mesh  intervals  in  each  direction  may  be  prescribed. 
Since  the  analysis  becomes  increasingly  unstable  as  the  mesh  is  refined,  it  is 
necessary  to  maintain  as  coarse  a  mesh  as  possible.  For  a  spherical  nose,  a 
5x5  T  -y  mesh  yielded  very  good  results.  In  general,  even  for  the  more  complex 
nose  geometries,  an  8  x  8  T  -y  mesh  should  be  a  rough  upper  limit  of  mesh  density. 
The  program  cannot  accommodate  a  mesh  more  refined  than  15  x  15. 

Columns  3-10  indicate  the  T  mesh  interval  adjacent  to  the  bow  shock 
(  T  =  1  ).,  and  cols.  43-50  refer  to  the  y  mesh  increment  adjacent  to  the  axis  of 
symmetry.  Refer  to  Fig.  II  and  the  section  on  sample  inputs  and  outputs  for  further 
details. 

Cards  C  +  1  and  C+  2  prescribe  a  new  set  of  "iteration"  points  which  are 
to  be  satisfied  by  the  perturbation  technique  at  the  conclusion  of  the  second  series 
of  sweeps  through  the  transonic  and  subsonic  flow'  regions.  These  points  need  not 
be  identical  to  those  used  initially;  the  number  of  points  may  also  differ. 
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Thus  the  number  of  input  cards  necessary  for  a  given  run  varies  from 
roughly  10  to  an  upper  limit  of  24;  in  most  cases  the  number  will  not  exceed  14. 

The  option  of  using  varying  mesh  intervals  was  included  so  that  a  relatively 
dense  mesh  may  be  prescribed  in  those  portions  of  the  subsonic  region  where  the 
velocity  gradients  are  high,  without  impairing  the  stability  (and  accuracy)  of  the 
solution.  A  recommended  mesh  for  a  spherical  nose  is: 

5  intervals  @  Ay  =  0.2 

5  intervals  @  =  0.2 

For  those  nose  geometries  which  exhibit  a  very  rapid  change  of  curvature 
in  the  neighborhood  of  the  sonic  point  (e.  g.  the  Apollo  configuration)  it  will  be 
necessary  to  utilize  3  mesh  intervals  in  each  direction  to  obtain  accurate  results. 

Note  that  the  last  "body  card"  (card  A)  should  be  blank.  Thus,  these 
may  be  up  to  12  body  profiles,  followed  by  a  blank  card,  followed  by  the  nose  geometry, 
etc.  The  value  of  Xj.  for  all  body  profiles  must  be  non-zero.  The  value  of  x.Q  should 
be  positive  (  a  0  )  on  card  A  +  1. 

TT-xpsri sn.c s  li3,s  ind.ic3.tsd.  th.2tt  no  moi* s  th.3n  3  1 1  it s 2* 2.ti on1 1  points  (including 
the  stagnation  point)  should  be  prescribed  for  the  first  series  of  sweeps  (cards  A+l, 
A+2,  and  A+3  ).  If  the  nose  is  spherical,  these  same  three  points  should  be  prescribed 
on  cards  C+l  and  C+2.  If  the  nose  is  not  spherical,  a  total  of  four  points  may  be 
prescribed  on  cards  C+l  and  C+2.  While  it  is  possible  to  prescribe  more  than  four 
points,  the  equations  on  page  32  tend  to  become  singular,  and  the  resulting  solution 
is  poor  -  thus  the  practical  upper  limit  is  four  iteration  points.  (The  example  listed 
on  page  S2  is  unfortunate  since  it  represents  an  early  test  case  which  conflicts  with 
the  above  rule.) 
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To  locate  these  iteration  points  properly  along  the  body  profile,  the 
following  approximate  criteria  are  recommended: 

If  we  define 

S  as  the  distance  from  the  axis  to  the  assumed  sonic  point,  along  the  profile 
of  the  body,  then  for  3  iteration  points, 


For  4  iteration  points 

S1  '  ° 

S  =  l/2  S 

S4  ■  S' 

Thus  for  a  spherical  ncse  of  unit  radius,  S  =  '4 

51  =  0,  Xj  =  yj  =  0  ( stagnation  point) 

52  =  W/6,  =  0,134  ,  yx  =  0.5 

53  =  ^4,  xl  =  0,293,  yl  =  0.  707,  for  3  points. 
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DESCRIPTION  OF  IBM  709/90/94 


COMPUTER  PROGRAMS  AND  ANALYSIS  FOR  FLOW  FIELDS 
ABOUT  BODIES  OF  REVOLUTION  IN  HYPERSONIC  FLIGHT 

By 

E.  Lieberman 


INTRODUCTION 


This  report  embodies  all  the  necessary  information  required*  for  the 
use  of  the  subject  programs.  It  is  assumed  that  the  user  has  a  general  background 
in  the  physics  described  herein,  and  it  is  recommended  that  the  analysis  be 
carefully  studied  before  any  attempt  is  made  to  utilize  these  programs.  While  a 
knowledge  of  computers  is  not  essential,  it  is  suggested  that  aid  should  be  sought 
where  necessary  to  ensure  that  the  input  data  cards  are  correctly  punched*  and 
that  the  operating  procedures  are  properly  executed. 

The  large  section  on  Flow  Diagrams  is  included  for  the  purpose  of 
completeness  and  is  intended  to  serve  as  a  guide,  only  for  those  expert  program¬ 
mers  who  have  the  coding  available,  and  who  wish  to  study  the  coding  in  some 


detail. 
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Description  of  Programs 

The  programs  described  in  this  report  make  it  possible  to  generate 
the  entire  flow  field  about  axisymmetric  bodies  of  revolution  in  hypersonic  flight. 
These  "second  generation"  programs  (see  Ref.  1)  can  accommodate  virtually  any 
geometrical  configuration,  including  those  having  very,  blunt  noses  and  after-bodies 
with  expansion  and  reentrant  corners.  The  gas  considered  is  in  chemical 
equilibrium  and  includes  the  effects  of  gas  dissociation  and  vibrational  excitation. 

These  programs  have  been  coded  for  the  IBM  709/7090/7094  large 
scale  digital  computers  having  a  32  K  (32,768  words)  core  capacity,  and  a 
minimum  of  4  tape  units  on  each  of  2  Data  Synchronizer  channels.  An  on-line 
card  punch  is  not  necessary.  Due  to  the  logical  nature  of  the  analysis  coupled 
with  storage  problems,  it  was  deemed  necessary  to  code  the  major  portion  of 
the  program  in  the  FAP  language.  Every  effort  was  made  to  minimize  the  cost 
of  using  these  programs;  the  entire  flow  field  can  be  generated  in  as  little  as 
5  minutes  on  the  IBM  7090  computer.  However,  for  some  configurations,  depend¬ 
ing  upon  the  free  stream  conditions,  bluntness  of  nose  region,  and  length  of 
after-body,  this  machine  time  could  increase  to  as  much  as  20  minutes;  average 
running  time  is  approximately  10  minutes. 

To  attain  the  utmost  in  efficiency  and  flexibility,  a  total  of  3  program 
decks  were  formed.  They  are: 

1.  TRANSONIC- SUBSONIC  program 

2.  SUPERSONIC  program 

3.  FLOW  FIELD  program 
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1.  TRANSONIC- SUBSONIC  Program 

This  program  accepts  as  input  data,  the  free  stream  conditions, 
coordinates  of  points  which  prescribe  the  geometry  of  the  nose  region  of  the  body, 

/ 

and  a  description  of  the  mesh  used  in  the  numerical  procedure  of  computing  the 
flow  properties  in  the  elliptic  region.  These  subsonic  flow  properties,  together 
with  the  flow  properties  at  points  along  a  second  family  characteristic  (reference) 
line  which  is  in  the  transonic  flow  region  just  downstream  of  the  sonic  line, 
constitute  the  output  of  this  program.  This  output  is  written  on  two  magnetic 
tapes:  one  is  a  BCD  tape  which  is  listed  on  peripheral  equipment,  and  the  other 
is  a  binary  tape  which  may  be  retained  for  subsequent  use  by  the  SUPERSONIC 
program.  The  procedure  followed  by  the  program  is  outlined  below: 

1.  By  satisfying  continuity  conditions  in  the  nose  region  of  the  body, 

the  approximate  geometry  of  the  detached  (bow)  shock  is  obtained  and  expressed 

2  4 

in  the  form,  x=a  +  a  y  +  a  y  .  For  some  nose  configurations,  this  form 
O  u  t: 

of  the  "basic"  shock  polynomial  will  be  inadequate,  resulting  in  envelopes  in  the 

8 

transonic  flow  region.  The  term  a  y  is  then  added  parametrically  to  overcome 

o 

this  situation. 

2.  The  program  then  marches  up  the  bow  shock  and  tonstructs  a 
mesh  of  characteristic  lines  row  by  row,  spanning  the  area  between  bow  shock 
and  body,  in  that  (transonic) region  of  the  flow  field  just  downstream  of  the  sonic 
line  (fig.  IV).  The  reference  line  is  located  and  transformed  into  the  j  -y  plane. 
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3.  With  the  geometry  of  the  bow  shock,  and  the  transformed 
properties  along  the  reference  line,  available,  the  calculations  proceed  by  an 
inverse  method,  to  determine  the  properties  within  the  subsonic  region.  The 
calculated  profile  of  the  nose  between  the  axis  of  symmetry  and  the  point  on  the 
body  common  to  the  reference  line,  is  Compared  With  the  given  profile*  Deviations 
in  the  x- direction,  between  the  two  profiles,  are  calculated  and  stored. 

4.  One  by  one,  the  coefficients  of  the  ''basic"  shock  polynominal 
obtained  in  step  1,  are  perturbed.  After  each  perturbation,  the  program  sweeps 
through  the  transonic  and  subsonic  regions,  as  described  in  steps  2  and  3, 
calculating  and  storing  the  resulting  deviations.  There  are  as  many  coefficients 
perturbed  and  resulting  sweeps,  as  there  are  prescribed  "iteration"  points  in 
the  nose  region. 

5. -  At  the  conclusion  of  this  procedure,  it  is  possible  to  solve  for 
the  corrections  to  be  applied  to  the  coefficients  of  the  "basic"  shock  polynomial 
to  form  a  "final"  polynomial,  such  that  agreement  is  attained  between  the 
calculated  and  the  prescribed  nose  profiles .  Using  this  "final"  polynomial  to 
describe  the  geometry  of  the  bow  shock,  the  program  sweeps  once  more  through 
the  transonic  and  the  subsonic  regions,  calculating  the  final  properties  along 
the  reference  line,  and  within  the  subsonic  region, 

6. -  The  free  stream  conditions,  final  reference  line,  and  body  profile 
equations,  are  Written  on  binary  tape  for  subsequent  use  if  desired. 
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2.  SUPERSONIC  Program 

This  program  accepts  as  input  data,  the  free  stream  conditions, 
properties  along  a  reference  line,  and  equations  which  prescribe  the  geometry 
of  the  body  profile,  from  the  refer ence  line,  downstream.  The  flow  properties 
along  the  bow  shock,  along  the  body,  and  within  the  flow  field  at  points  on  the 
mesh  of  characteristics,  constitute  the  output  of  the  program.  An  "up-dated" 
binary  tape  is  also  written.  The  procedure  followed  by  the  program  is  outlined 
below. 

1.  March  up  first  family  characteristic  lines  (slope  =  tan  (0  +  ^  ) 
from  the  reference  line  to  the  bow  shock,  starting  near  the  intersection  of  the 
reference  line  and  bow  shock,  and  gradually  extending  the  domain  until  region  A 
is  completed  (see  fig.  I). 

2.  March  down  second  family  characteristic  lines  (slope  =  tan (0- (j)  ) 
from  the  final  first  family  line  in  region  A  to  the  body,  starting  near  the  inter¬ 
section  of  this  first  family  line  and  the  body.  Eventually,  region  B  of  the  flow 
field  is  completed. 

3.  The  remainder  of  the  flow  field  is  denoted  as  region  C.  This 
region  is  constructed  by  first  computing  a  point  on  the  bow  shock  by  utilizing  the 
first  two  points  on  the  previous  second  family  line,  and  then  marching  down  a 
second  family  line  from  this  shock  point,  to  the  body. 

4.  The  logic  used  in  regions  B  and  C,  permit  the  program  to  detect 
discontinuities  in  the  body  profile.  Thus  an  expansion  corner  -  which  produces  a 
Prandtl  -  Meyer  fan  -  or  a  reentrant  corner  -  which  is  the  origin  of  a  secondary 
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shock  -  may  be  prescribed  at  any  point  downstream  of  the  reference  line. 

5.  The  program  can  accommodate  any  combination  of  expansion 
corners  and  reentrant  corners  subject  to  the  restriction  that  no  second  family 
characteristic  line  can  cross  more  than  two  secondary  shocks  prior  to  intersecting 
the  body.  The  program  recognizes  that  a  secondary  shock  may  "die"  within  the 
interior  of  the  flow  field,  and  the  logic  adjusts  itself  accordingly.  Intersection 

of  two  secondary  shocks  within  the  flow  field,  or  of  a  secondary  shock  and  the 
bow  shock,  will  be  correctly  calculated  and  the  mesh  adjusted,  only  if  one  of  the 
two  intersecting  shock  is  very  much  "weaker"  than  the  other,  at  the  point  of 
intersection.  If  the  two  shocks  are  of  substantially  the  same  strength,  then  the 
resulting  shear  layer  creates  a  discontinuity  in  the  mesh  downstream  of  the 
point  of  intersection,  and  this  condition  lies  outside  the  scope  of  this  program. 
Another  condition  which  could  arise  is  the  case  of  a  strong  secondary  shock  causing 
the  flow  downstream  to  become  subsonic,  in  which  case,  the  method  of  character¬ 
istics  does,  not  apply. 

Another  restrictive  condition  is  a  situation  where  the  slope  of  a 
second  family  line  is  essentially  the  same  as  the  slope  of  a  portion  of  the  body 
profile.  This  would  create  a  large  gap  in  the  mesh  of  characteristics  which 
could  cause  the  program  to  abort.  While  there  are  logical  mechanisms  built 
into  the  coding  to  control  the  mesh  to  some  extent,  an  intrinsic  condition  such 
as  that  just  described,  can  not  be  remedied.' 

6.  Bodies  with  conical  noses  fall  within  the  scope  of  this  program; 
of  course,  only  the  supersonic  flow  field  exists  for  these  configurations.  A 
subroutine  which  is  included  as  part  of  the  SUPERSONIC  program,  calculates 
properties  along  a  horizontal  reference  line,  and  thus  serves  as  a  prologue  to  the 
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major  portion  of  the  program. 

3.  FLOW  FIELD  Program 

This  program  is  a  CHAIN  job  consisting  of  two  LINKS.  Essentially, 
the  first  link  is  the  TRANSONIC -SUB SONIC  program,  while  the  second  link  is 
the  SUPERSONIC  program.  The  automatic  linkage  of  the  two  parts  is  accom¬ 
plished  by  utilizing  a  feature  of  the  IB  FORTRAN  Monitor  System  and  the  transfer 
of  common  data  is  done  via  the  binary  tape  previously  described. 

The  input  to  this  program  is  identical  to  the'  input  required  by  the 
TRANSONIC- SUBSONIC  program,  and  the  output  is  the  sum  total  of  the  two 
components  of  the  FLOW  FIELD  program. 

The  above  description  merely  presents  the  overall  scope  of  application 
of  the  subject  programs.  All  necessary,  details  are  presented  in  the  following 
sections. 

While  these  programs  have  been  tested  for  a  wide  range  of  flight 
conditions  and  body  configurations,  it  is  certain  that  not  all  possible  conditions 
and  permutations  of  branching  have  been  encountered.  In  addition,  due  to  the 
limitations  of  time  and  cost,  some  of  the  more  "marginal"  features  of  the 
SUPERSONIC  program  (LINK  2),  such  as  shock  intersections  (which  in  general 
lies  outside  the  scope  of  this  program  as  indicated  previously),  and  certain 
variations  of  "dying"  secondary  shocks,  have  not  been  as  vigorously  tested  as 
the  major  components  of  the  program. 


CALCULATION  OF  COEFFICIENTS  OF  "BASIC"  POLYNOMIAL  TO 


APPROXIMATE  DETACHED  SHOCK 


A  rough  guess  of  the  detachment  distance,  5q  ,  is  found  as  follows: 

The  radius  of  curvature  of  the  body  at  the  "sonic"  point  (where  slope 

* 

is  unity)  is  calculated  numerically  and  denoted  by  , 


DISC 
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Next,  cross  a  normal  shock  to  obtain  the  density,  6  ,  and  y ,  behind  it  and 

"  n 

calculate. 


Sphere: 


^o_ 

h 
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-1 


and 


Disc: 
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T  *  >*03 


00 
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*  .  o  * 

where  h  =  h  +  R  (  1  -  cos  45  ).  Note  that  h  =  y 

s  s  s 

R  * 

Then,  for  the  given  body,  s/  h  may  be  calculated  and  the  corresponding  value  of 
o  /  h  may  be  found  by  .linearly  interpolating,  between  those  values  calculated  for 
the  disc  and  for  the  spher^,  as  indicated  in  the  sketch.  Thus,  multiplication  by 
h  yields  the  initial  estimate  of  the  detachment  distance,  ®o,  for  this  body.  Implied 
in  this  approach,  is  the  assumption  that  the  point  (  x  ,  y  )  closely  approximates 
the  sonic  point  on  the  bodjr. 

Next,  we  improve  the  approximatidn  of  $>,  iteratively,  as  follows: 

Compute 
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where  is  the  radius  of  the  body  at  the  axis  of  symmetry,  calculated  numerically. 
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Then 


dto\  =  _  C_ 
~dsb  B 


(3)  and  we  may  calculate  the  radius  of  curvature  of  the 


detached  shock  at  the  axis,  R 
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Solve  the  following  two  equations  for  the  coordinates  of  points  s: 
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X  =  x  +  y  -  y 


2  $  * 
y  +  2  R  y  =  2  R  (x+y-x+6) 
8  sl^  S  shQ  o  o 

Note  that  the  first  two  terms  of  the  detached  shock  polynominal  are 

ao  =V«,>and  S  =(2R-u  )_1 
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To  compute  the  mass  flow  across  the  shock  layer,  6  ,  it  is  necessary  to  compute 

the  pressure,  ,  and  velocity,  v^  ,  behind  the  shock  at  point  s.  In  addition 
these  properties  (  ,  v^'  )  must  also  be  computed  at  the  sonic  point  on  the  body. 
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It  is  necessary  to  compute  the  rate  of  change  of  mass  flow  in  the  direction  normal 
to  the  body  at  the  sonic  point,  and  then  assume  that  the  variation  of  mass. flow  is 
closely  approximated  by  a  parabolic  across  the  shock  6^  .  Then 
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where  n  is  the  normal  to  the  body  at  point  r 
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is  the  radius  of  curvature  of  the  body  at  point  r  H  Rg 
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i  W 
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\  y 
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■  y  n 


wr 
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J 
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a  "  ^  v  - 1  ^  ^  y 
\  ^  Tb  7 


where 
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I 
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I 
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r 

c 
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Fronj  Bernoulli's  equation  we  solve  for 


Then  assume  the  variation  along  6 


L 

[ 

[ 

I 

[ 

1 

i 
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Then  A  =  P  v 

2  wo  w„  b  b 
s  s 


*  . 


Integrating  across  this  shock  layer,  the  mass  flow,  M  ,  is  calculated. 


M 


* 

Pb  Vb  , 
2  cos  45° 


[v  +6 

Lr  r 


cos  45 


-  y_ 


+  A  6 
l  r 


?  i  i\  y 

Compare  M  =  „.!L.  ~  1  _ SL 

by  +  l  2 


V  0  cos 

+  -£ - 


45? 


y  6  cos  45°  I 


K-. 


'y  o 

Az6r-^+- 
rL  3 


(9) 


If  M 


HillI 


y  +  1 


y  2 

s 


then  6=5  + 

o  o  — 


A  6 


where  A6q  is  a  prescribed  small  increment.  Return  to  eq.  (1)  with  this  new 
value  of  6q  ,  compute  new  coefficients  aQ  and  a2  and  check  M  again  using  (9). 
Continue  this  iteration,  adjusting  6q  as  required,  until  continuity  at  the  shock 
layer,  6r  ,  is  satisfied.  Then  test 


As  < 

— -  >0.5  If  <  0.  5,  the  basic  shock  polynomial  must  be  refined  by 


adding  an  y  term.  The  other  condition  lies  outside  the 
scope  of  this  program. 

The  coefficients  aQ  and  a2  are  held  fixed  while  a^  is  determined  by 
satisfying  continuity  across  the  shock  layer,  6,  .  Using  y  across  the  normal  shock, 


initially  set 

u  =  3  v  -  e 


-  J111 

2  V  y  -  1 


.tan 


-1 


/  r 

1  "1 
L 

|2  ( 

r  ~  1 

nc 

y; 

]  , 

-  i 

|  2  ‘ 

2 

/■« c\ 

r 

-l 

r 

/ 

n 


is  approximated  by  2  sin  0  . 


where 


/ 


Then  iterate  on  jac  ,  solving  for  Pt  and  vc  in  the  process 
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r  rJ  » 

Xll 

y 

=  1/2' 

K  J 

1  +  cos  2' 


%4l  14 


■0-  "He 


(10) 


n_  €_  3 it_ 

since  Gr  4  pr  =  4  +  2  =  4 


£  X  n 

Pc*4  » 


1  ly 


c 

m 


2  =  c  -  2y  nc 


vc  =c 


Y  +  l  K~ 


(ID 


and  jic  =  sin 


-1 


Jyi 

rJ  ^ 

nc 

141 

y  y+i 

^c  - 

Vc 

>  =  tan 


-1 


_Y 

xzl 

y.4  1 

rJ 

n  . 

*!h. 

"4^ 

O 

N. 

1 

Y-l 

Y  41 1 

1 

1  i 

[cj 

1  /  2 

Compare  this  value  pf  p,c  with  that  used  in  (10) ;  return  to  (  lQ  with  this  value 

(n+1)  (ri)  I 

if  no  agreement  is  reached  arid  repeat  until  |id  -  p.c  |  10 


-4 


It  is  now  possible  to  calculate  thd  value  h4  of  the  shock  polynomials  . 

,  4  , 

x=  aQ+  a2y  +  a^y  ,  by  satisfying  continuity  across  the  shock  layer,  6C.  At  point  d 

4 


xd  -  ao  -  a2yd  -  a4  y^  =  0 

Setting  5q=  6r  initially,  then  with  xc  ,  yc  ,  0C  ,  given  (inputs), 

Yd  =  Yc  +  6c  co8  ®c 

xd  =  xc  "  6c  ,ln  ®c 
1 


(12) 


a4  =  - 


dx 


yd 


a2  yd  +  »o  -  xd 


U3) 


dy  *  tan  «  d  =  2  a2yd  +  4  a4yd 
'd 
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Compute  Pw  and  vw  from  (  5)  and  then  the  mass  flow  across  the  shock 
layer,  6C,  may  be  computed.  It  is  necessary  to  compute  the  rate  of  change  of 
mass  flow  in  the  direction  normal  to  the  body  at  point  c,  and  then  assume  that 
the  variation  of  mass  flow  is  closely  approximated  by  a  cubic  law  across  the 
shock  layer,  6C.  Then 


— )  N  +  1  \  Pcvc2 

1  ^n  c  lv”l  I  Rc 


where  n  is  the  normal  to  the  body  at  point  c. 

is  the  radius  of  curvature  of  the  body  at  point  c,  given  (input). 
Pc  ,  vc  are  calculated  from  ( 11)  . 


TTVorri  e V*r» r V  rplatinnR 


/ 


(15) 


From  Bernoulli 


s  equation,  we  solve  for 


2y 

—  1 

1 

y+lj 

^T. 

'/Vi 

an 

8n  c 


and 


9(  Pv) 

—  TD 

I  By 

+  VC 

r 

j  aJ 

ap 

dn 

-  Pc 

(  an 

(  8n 

- 

c 

> 
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(16) 


Then  assuming  th<4  variation  along  6C, 


/>J 


Pv  =  pc  vc  +  Ai 


n 


n 


j+  A3  I  g  J-  f  solve  for  A3 
A3  "  Pwd  VW(J  “  pcvc"  A1 


Integrating  across  the  shock  layer, 
6 


M 


r  ^ 

=  I  Pv  (  yc  +  n  cos  0C)  dn 


or 


M 


M= 


J 


pcvc  +  Ai 


6  + A: 


n 


yc+  6C 


n 


CO'S  Qr 


6cdIT 


or 


L 

J  L 

J 

/A/ 

“ 

- 

“  ,  " 

pcvc  ( 

y c+  cos  Qc 

2  2 

^  +  A1  5C 

Y-c  6C  cos  Oc 

+a36c 

Yc  t  6e  mg  Oc 

2cos  ec  | 

2  +  3 

4  5  • 

1 

- 

.  J 

(  17) 


Compare  M  = 


r-1  1  yd 


y+  1  j  2 
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If  M  ^  ^77 
•y  +  1 


1± 

2 


then  +  A  5C 


where  A  is  a  prescribed,  small  increment.  Return  to  (  13)  ,  (  5) ,  then  (  14) 
through  (  17)  ,  and  iterate  on  6C  until  continuity  is  satisfied.  After  convergence, 
the  value  of  a^  from  (10)  is  held  fixed,  and  from  a  system  analogous  to  (  12)  , 

Ys=  Yr+  ^r  cos  45  ,  and  since  tan  0r  -  1 » 

(18) 

4  2 

a4y8  +  ^  Ye  +Ys-yr+ao-xr  =  0 


Solve  (  lg)  iteratively  for  ys,  adjusting  6r  until  (18)  is  satisfied.  Then, 

3 

tan  w  =  2  a2  y  s  +  4  a  .y  s 
s 

Compute  (5)  through  (  9)  ,  testing  continuity  and  adjusting  bQ  if  not  satisfied; 
execute  eqs.  (1)  through  (  3a)  .  Then  with  aD=  xD  -  6Q  and  a2=  (  2  RshQ)  >  return 
to  (  18)  for  6r  and  ,  Continue  this  iteration  on  6Q  until  continuity  is  satisfied. 

When  continuity  is  satisfied  across  the  shock  layer,  6r ,  hold  aQ  and  a2  fixed, 
and  return  to  (12)  to  iterate  again  on  6C,  revising  a^  in  the  process.  And  so  the 
double  iteration  to  determine  those  values  of  aQ,  az,  and  ag  which  define  a  basic 
shock  polynomial  continues  until 


M, 


II 

LTLlIli] 

Y+  1/  2 

(  1  +  0.025)  , 


i.e.  continuity  is  satisfied  within  23%  across  both  shock  layer  considered, 
simultaneously. 

For  certain  nose  configurations  which  differ  significantly  from  spheres,  it  is 
necessary  to  prescribe  the  downstream  portion  of  the  shock  geometry  with  great 
precision  to  successfully  generate  a  reference  line  and  avoid  the  possibility  of  an 
envelope  within  the  flow  field.  In  such  cases  it  is  necessary  to  add  to  the  basic 
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shock  polynomial,  the  term  agy  .  Since  the  coefficient,  ag,  is  extremely 
small  (  and  negative)  ,  this  term  has  virtually  no  effect  upon  the  mass  flow  cal¬ 
culations,  and  is  merely  used,  parametrically,  as  a  device  to  successfully 
execute  the  transonic  region  of  the  flow  field. 

Note  that  the  preceeding  analysis  is  based  upon  the  assumptions 
of  linear  mass  flow  variation  between  shock  and  body  at  the  axis,  a  parabolic 
variation  across  the  shock  layer,  r  s,  and  a  cubic  variation  at  the  downstream 
section,  cd.  These  assumptions  are  consistant  with  the  physics:  as  one 
proceeds  around  the  "shoulder",  in  the  downstream  direction  this  variation 
assumes  still  higher  orders. 
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TRANSONIC  ANALYSIS 


After  establishing  the  coefficients  of  the  polynomial  which  prescribes  the 
shape  of  the  detached  shock  in  the  nose  region,  the  program  "marches  up"  this 
shock  from  the  axis  until  it  finds  the  first  supersonic  point.  It  then  continues 
its  march  downstream,  generating  properties  behind  the  shock  at  S3  additional 
points,  constructing  a  mass  flow-entropy  table  in  the  process. 

With  this  information,  a  mesh  of  characteristic?  is  constructed,  extending 
from  shock  to  body  as  indicated  in  Fig.  IV.  The  reference  line  is  that  second 
family  characteristic  line  which  extends  from  the  shock  and  terminates  at  the 
supersonic  body  point  which  is  farthest  upstream  (closest  to  the  axis)  . 

At  the  conclusion  of  the  Transonic  analysis,  the  program  executes  a 
routine  which  transforms  the  properties  along  the  reference  line  from  the  physical 
plane  to  the  T  -y  plane  which  is  used  by  the  Subsonic  analysis.  Some  small 
inconsistencies  arise  due  to  the  fact  that  y  varies  from  point  to  point  in  the 
Transonic  region,  while  it  is  assumed  constant  in  the  Subsonic  analysis.  These 
inconsistencies  are  minimized  by  assigning  as  the  y  used  in  the  Suhsonic  region, 
the  average  of  the  values  at  the  body  point  and  at  the  axis  behind  the  detached 
shock  (  across  a  normal  shock)  .  These  inconsistencies  are  adjusted  along  the 
length  of  the  reference  line,  and  it  has  been  found  that  they  have  no  noticeable 


effect  on  the  results. 
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TRANSONIC  ANALYSIS 


The  properties  behind  the  detached  shock  at  54  points  are  calculated 
in  the  same  manner  as  those  in  the  Supersonic  program  (  see  page  37).  The 
only  difference  is  that  the  value  of  e  at  every  point  is  calculated  from  the 
derivative  of  the  shock  polynomial,  i"  (y),  rather  than  by  computing  A  €>  as  Is 
necessary  in  the  supersonic  program. 

INTERIOR  POINT 
The  properties  of  points 
A,  B,  and  coordinates  of  DB 
are  known;  the  properties  of  point 
C,  and  the  qoordinates  of  point  DC  are  required. 

Points  DB  and  DC  are  defined  as  those  shock  points  having  the  same  entropy  as 
points  B  and  C,  respectively. 


x  = 
c 


yB  ~yA'  XB  tan(  }B  +  XA  taii(  >A 

tan  {  0+  u  )  .  -  tan  (  0-  fi  )_ 

A  Jd 


yc  =  yA+  (XC-XA}  tan(e  +M  }A 


DC 


[ 


y 2 

DB 


%  WI 


(  y l  -  y4  ) sin  m 


B 


p  w 

00  CO 


sin  (  0  -  fi  ) 


XL 


B  J 


x  =  F  (  y  ) 
DC  DC 
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If  y  jjq  <  0»  then  point  C  lie*  within  the  domain  of  the  body  -  see  analysis 
of  Body  Point,  below. 


extract 


\R/ 


Compute  the  mass  flow  at  point  DC  and  enter  the  entropy  table  to 

'  Si 


iR; 


Then 


DC 


S 

iTTl 


DC 


The  remaining  properties  are  calculated  in  the  same  manner  as  an 
interior  point  for  the  Supersonic  program  (page  43  ). 

BODY  POINT 

Since  the  test  v  <0 
y  DC 

identifies  point  C  as  being 
within  the  domain  of  the  body, 
it  is  necessary  to  literate  to  locate  a  point  on  the  surface  of  the  body,  which 
corresponds  to  y  qq  =  0. 

The  program  divides  the  line  BA  into  4  0  segments,  establishes  the 

properties  of  points  1,2,  3, . by  interpolation,  and  successively  calculates 

new  points  "C",  using  point  B  and  points  1,  2,  3,  .... .  until  y  _  <  0  once 
again.  By  interpolating,  the  properties  of  point  C  ,  on  the  surface  of  the 
body  are  determined  within  very  close  tolerances. 


SUBSONIC  ANALYSIS 
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At  the  completion  of  the  transonic  region,  the  reference  line  is  defined, 
and  the  program  then  transforms  the  properties  of  all  points  on  the  reference 
line,  from  the  physical  x-y  plane  to  the  T -y  plane.  The  properties  along  the 

reference  line  required  by  the  Subsonic  analysis  are  U,  V,T  and  y, 

W  W_ 

U  =  yy^cos  9  ,  V  =  yy^  sin  6  ,  y  =  y 

and  T  is  found  by  integrating  along  the  reference  line  from  the  shock  (T«  1) 
towards  the  body  according  to  the  following  formula.  Refer  to  Fig.  II. 


^.n'^ns 


-n+2ynq--(y  n-  Vue)  +[(R-q)n+(R-U)nB 


(  y  n";  yns>  -(xn“xns) 


(R.V)n+(R.y 


JfVn+ynsJ  , 

[  2  +  (yn-YnS) 


where  6  =  (  ::  p-  5 


y  +  1 


n 


CD 


y  =  i/2  (yb+  yn) 


yb  is  the  value  at  the  body  on  the  reference  line 

y  is  the  value  behind  a  normal  shock  (  detached  shock  at  the  axis) 
n 


Since  .^body  not  equal  zero  identically,  the  "noise"is  distributed 

uniformly  along  the  reference  line  as  follows: 

T  =  Tn'-^b 
1-^b 

where  Tn*  is  the  value  calculated  by  the  integration  scheme,  above. 

T"b  is  the  value  calculated  at  the  body  (  / 0) 


)s  the  adjusted  value  at  each  point,  n,  on  the  reference  line. 

Tb  is  assumed  very  small  compared  to  unity  ( «0.  1)  ,  and  test  runs  have 
indicated  that  this  assumption  is  correct. 


hs 
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SUBSONIC  ANALYSIS 


At  each  interior  mesh  line  of  constant  T,  we  set  up  and  solve  a  system  of 
linear  equations  in  the  unknown  velocity  components,  U  and  V,  Since  the  analysis 
is  non-linear,  the  coefficients  of  the  system  being  functions  of  the  unknowns,  it  is 
necessary  to  approximate  the  unknowns  first,  using  differencing  techniques, 
compute  these  coefficients,  and  then  solve  the  system  to  obtain. the  variation  of 
U  and  V  ( as  well  as  the  other  desired  properties)  along  the  mesh  line,  from  the 
axis  of  symmetry  to  the  reference  line.  To  refine  and  improve  the  stability  of  the 
solution^  the  procedure  is  repeated  once  more,  utilising  the  more  accurate  values 
of  U  and  V  resulting  from  the  first  solution.  Notably,  the  derivatives  with  respect 
to  T  are  significantly  refined.  Refer  to  Fig.  II. 

Yn  “Yns  v. 

tan  ‘t’ns  =  | 


dU  ■ 
dT  , 


dV 

dT 


ns 


9U  9U_ 
9T  +  9y 


tan  <j) 


ns 


9V  ,  9V  , 

8T  +  57" 


ns 


ns 


Fbr  a  first  approximation  of  the  unknowns  Un  and  Vn., 
un  =:  uns  +  (5TL(Tn-T„B) 


ns 


Vn  “  Vns  +  ("dT^s  (Tn  ”  Tns) 


(  1) 


Substitute  the  above  Un  and  Vn  into  the  following  Equation  (  2) 
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where  f  sf(  Ty2)  =  P*R*^Y)  6  =  Y  ~  1 
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c  = 


(Ym‘V  K 


00  1  CO  J 


+  1 


1 

— -  + 

1 

1 

a2  (  y  + 1 )  +  2,3  ,  1 

V  M  2 

1 00  00 

( Y+  1 )  (  1  +  g) 

T 

(l  +  g)  ( Y  +  if  d+gf 

*1  - 

1/2 


(V+  1) 


+  a  (  1+  g) 


(7+T)2  '  (P^TT  + 
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8T  6  vir  8T  ' 
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( Y  +  1 ) 


+  P  )  -a2  (1+  g) 


1  (,7^r+2p)(l+g)  +  (/(l+g)^ 


(v  + 1) 2  M\+i) 


+  (4) Y 

R  (l  +  g)2  9t 


<-• 


1 


(Y  +  1) 


+  f»-  1 


(y+D  (v+i)^(i+g^ 

T  7 


-  (—  +2(3) 

- 

(l+g) 


l 

T7U 


+  (Y  +  l.)Z(l+g)2 


-£ 


m 


mi 


m=  2,  4 


CmT  •  -2  '  m  and  M  =V^  —  C  T 
y  m 


<T- 1) 


rm 


where  g  =  g(  T>  y)  >  consisting  of  the  terms  of  the  polynomial,  Fy  ,  with  the 
T  term  inserted  as  indicated. 

F(  y)  is  polynomial  defining  shock  curve 
Fy  =  dF/dy 


Y  +  1  v  Jvi  2 

•  OO  00 


(y-D  f(  ygo 


P  =  t'V+  1)  ['(Voo-1)  (Y-l)  7  Yoo^Jj 

\ 

^in  ~  ^in  "  ^in  tan  ^  ns 


D.  =  D.  -  B.  tan  4> 
in  in  in  ns 


1/2 


I 


i  =  1,2 


At  point  n,  two  equations  are  obtained  (Unip  -1) 
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ain  Un  +  cin(^i+  tfU(n-l))  +  binVn+  dinC^h+  ifV(n-^  +  ein=  0 


i  =  1,2;  in  the  unknowns  U  and  V 


where  a^n  =  2A^n 


bin  =2Bin 


_  ( ^"n  “Tns) 

Cin  =  Cin  (^n+ij^n-i)> 


din  =  Din  ( 


(  ^n  "^ns  ) 

^jn+  lj-y(n-lj) 


dU  dV 

ein  =  “  -A-in  ^ns+  ®in  ^ns)  “■ (-^-in .(  +  ®in  (<jt  )  ~^in  )  (^n'^ns) 


When  n  =  1,  then  Un_i  =  U  J  is  found: 


y=  o 


(1-6)  . 

U  [(c-U2) /(  1+6)  f(oj  2A  =  T 6 


and  3V/0T  1 

J  y=  o 


atj 

ay  0=  v°=  0 


At  the  axis, 


U 


z  6-1 

_  T  fi  I”  (C-U  )  "1  ->  c 

y=  °  Ln  +  sywJ 


iteratively 

T 

initially  sejt  U  =■  — ^ 


(3) 
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dU  f dU \  ,  dY  /dVi 

For  the  first  solution,  ^  =  Idf  I  *nd  ^  j 


ns. 


However,  for  the  second,  "refined"  solution,  4~  =1/2 

dT 


|  du| 

JdU\ 

dT 

1  dT  1 

’  'ns 

.  dV 

and  —  =1/2 

dr 


/  dV'1 

/  dV  \ 

.  /  dU 

and 

dV 

dT 

+  dT 

where  — 

dr  , 

ns  'n 

» 

n 

are  calculated 


from  the  results  of  Eqs.  (  5)  evaluated  at  the  $fid  of  the  first  solution. 

Equation  (3)  forms  a  system  of  2(p-l)  linear  algebraic  equations  which 
are  solved  simultaneously  yielding  the  values  of  ,  Vn(l<ns  p-1)  .  Using 
these  values  the  following  expressions  are  evaluated: 


/SU|  _  U(  n-f  1  )~U  (  n-1 ) 

-  V+l)  V(n-1) 

Sy)»  y(n+l)-y(„-l) 


>  1  s  n  s  p-1 


(4) 
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The  solution  of  the  2  (p-1)  linear  simultaneous  equations,  and 
Eq.  (4)  and  (5)  above,  yield  the  necessary  parameters  at  every  point,  n, 
n  =  1,2, . p-1 . 

3U  9  V 

Using  the  values  of  Un>  ,  (  gj"  )n  and  (  )n  just  calculated, 
return  to  Eqs.  (  2)  ,  (  3)  ,  (  4)  ,  and  (5) .  Then  the  following  is  calculated, 
and  the  program  "marches"  to  the  next  mesh  line. 

To  affect  the  transformation  from  T  to  the  physical  x: 

Within  the  flow  field, 


dx 


Tv  [<vf> 


dT  +(T  -  6  -R  . 


U)  dy] 


in  difference  form: 
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At  the  axis , 


(  yn-yns) 


(  T  n“T  ns)  + 


(  T n+  Tns  )  6 


dx*  - 


U 

2T 


dT 


8V 
9y 

In  difference  form: 


i  u  )  (_H_) 

f  1  ns  f  1  n 

.ay.  .  av. 
(z~)  +(  t-  ) 
8y  ns  &y  n 


(  T  n-T ns)  “  xns‘ 


(Un+Uns)  ^  Tn”  Tns) 

(tTt  )  (  \ 

n  nsM9yn,  +  9yn 


Knowing  the  values  of  x  along  the  shock,  the  values  of  x  along  line  1  may  be 
computed,  and  this  procedure  repeated  until  the  body  is  described. 
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I 

I 

I 

I 

I 

[ 

[ 


L 

IT 

*  •! 

1 

1 1 


I 

I 

I 


BEHIND  SHOCK  WAVE 


U  =  1  - 

r  . ..  j 

2a 

d-  (y  +  i)  +  ,  i 

(y  +  l)  ( 1+  Fy‘) 

(1+Fy2)  '  (  1+ Fy2)2  (y +T)2 

1/2 


2a 


V  =  F, 


■M  1)(1+  Fy‘)  ”Q 


+  F„ 


(y  + 1) 


+  2(3 


(  1  +  FyT  +  (  1  +  Fyz )  2  ( y  +  1 ) 


1 

ITT; 


1/2 


»U  ZFyFyy 


3y  (1+F/)2 


+P)- 


l _  .  I  y  +  i  (y  +  i)’2(  i+  fv2) 


(y  +  l) 


2  a 

- ;  +  26 

a3--  y  +  1  + 


1  1/2 


(!+FyM  (y+  l)2  (  1+  Fy2  )2 


6y  =  Fyy  (1'u)‘  Fy  ^  '■'h*r',  Fyy= 

9U  9V 

Solve  the  following  Z  sim,  equ.  for  -  and - 

9t  9t 


A  3U  _  av  „  9U  ^  3V  ^ 

Ai  9T  +  Bi  9T  +  ci  ay  +  Di  9y  +  ex  =  o 


9U  9V  ,  r  9U  A  .  n. 

3T  9T  3y  2 


where  Aj  ,  ,  D^,  Ej  (i=  1,2)  are  found  from  (  Z)  . 

These  properties  are  calculated  at  mesh  points  along  the  mesh  line,  t-  1, 
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ITERATION  ON  DETACHED  SHOCK  POLYNOMIAL 


After  the  "basic"  shock  polynomial,  which  is  of  the  form,  x=  F(  y)  =  aQ  + 
2  4g 

a2  y  +  a  y  +  a  y  ,  has  been  generated  (  see  page  8  )  ,  the  program 

sweeps  through  the  Transonic  and  Subsonic  regions,  calculating  as  an  end 

result,  the  profile  of  the  body  in  the  vicinity  of  the  nose.  The  profile 

calculated  will  not,  in. general,  agree  with  the  prescribed  body  geometry. 

The  deviations  in  the  x  -  direction  along  the  body,  at  the  prescribed 

interation  points  (see  section  on  Input  Format,  page  77)  ,  i,  are  calculated 

( o) 

as  6  .  The  superscript  (o),  identifies  these  6  as  being  caused  by 

1  i 

the  "basic"  shock  polynomial.  See  sketch  below. 


AXIS  or  SY MM. 
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If  a  total  of  J  points  are  prescribed  as  iteration  points  (  excluding  the 

stagnation  point  -  J  =  3  in  above  sketch)  ,  the  program  executes  J  additional 

sweeps  through  the  Transonic -Subsonic  regions.  For  each  sweep,  one  coefficient 

of  the  "basic"  shock  polynomial  is  perturbed  and  the  resulting  deviations  calculated. 

Thus  for  the  jth  sweep,  j=  1,2,....,!,  the  detached  shock  polynomial  is  of  the 

form  F  (  y)  =  a  +  a2  y  +  a  y  +  a  y  +  A  a2;  yZJ,  and  the  deviations  are 
^  4  8  ^ 

prescribed  as  fi/ ^  ,  i  =  Note  that  a  maybe  zero, 

i  8 

After  J  sweeps  have  been  completed,  a  system  of  simultaneous  equations 


may  be  written  in  the  unknowns,  C.,  as  follows: 

Cj  61(1)  +  C2  6l<2>  +  C3  6i(3)  + 

Ci  62(1  }  +  C2  fi2(  2)  +  C3  62(3  }  + 

Cj  6j  (l)+  C26j(2)  +  C36j(3)  + 


V 


+  Cj6l  (J)= 

+  Cj62<J^  -62(0) 

,  r  ,  (J).  .  (0) 

+  Cj  6J  ~  '  6J 


-  z  4 

F  ( y)  =  ao  +  a2  y  +  ^  y  +  . 


zJ 


where 


r 


The  final  polynomial  is  then  written 

4 

+  a3j  y 

i2j  =  a2j  +  C.  .  Aa2j  ;  A  a2j  is  the  jth  perturbation. 


and 


j  =  1,2,3, 


a  =  a  - 
o  o 


,  J 


«  <0)  +  f 

c 

Vj)  -s<0,1 

°  M 

j 

A  final  sweep  through  the  T ransonic -Subsonic  region  uses  F  (  y)  to  describe 
the  detached  shock.  The  points  along  the  reference  line  are  written  on  binary^tape 
B3  along  with  the  free  stream  conditions  and  the  body  equations,  and  the  properties 
within  the  Subsonic  region  are  printed  as  output. 
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SUPERSONIC  ANALYSIS 


The  flow  properties  at  points  along  that  reference  line  generated 
during  the  final  sweep  through  the  Transonic  region,  represent  the  initial 
conditions  for  calculation  of  the  Supersonic  flow  field.  In  addition,  the  free 
stream  conditions  and  the  geometry  of  the  after-body  are  required. 

With  this  information,  a  mesh  of  characteristic  lines  is  constructed, 
and  the  flow  properties  at  each  mesh  point  are  calculated.  To  conserve 
machine  time  (  and  paper  ),  the  flow  properties  at  every  other  interior  point 
in  region  C  (see  fig.  1)  is  written  on  magnetic  tape  for  subsequent  listing. 
However,  if  a  secondary  shock  is  encountered,  thereafter  all  points  are 


written. 
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Preliminary  Properties: 


W  =  M  a  =  M  J  y  •  R  ■  T  1 

OD  00  OO  00  00  00 


K,  =  2  Cp  T  +  W 

1  ^  OO  00  00 


where  R  =  1716,  Cp  =  6006,  y  -  1.4 

oo  oo 


T  = 
t 


K, 


2  Cp 


oo 


r  -i 
00 

y 

00 


.i 


solve  by  iteration  where  Tt  initially  = 


K] 

2  Cp 


CD 


0  /  T 


E  =  e 
oo  — 


oo 


9'/T 


roo  = 


00 


oo 


-1 


y  -1 
oo 


I 

=  5500 
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Point  on  Detached  Shock 

The  following  properties  are  known 

at  Points  A  &  B: 

g 

x,  y,  P'.a  -  .h'.p'.y.U.m/w 

In  addition,  c  and  r  are  known. 

B  B 


The  object  is  to  generate  the  same 
properties  at  point  C,  as  at  point  B. 

To  determine  the  change  in  curvature,  of  the  shock  between  points 
B  and  C,  certain  derivatives  must  be  evaluated  at  point  B: 


Calling  =  cot 

p  =  0.002498 
B 


€ 

B 

2.302585  p' 
B 


P' 

=  2116.4  e  B 


u  =  W 
oo 


r  i  - (i  - 


oo 


B 


(  1  +  F4  ) 
y 


F 

y 


oo 


(1'pb  ’ 


ip 

d  € 


P  W  {  1  + 
co  "oo  ' 


Fy>'1 


PB 


2  ( 


F 

y 


1  d_p 

P  W  d  < 

00  oo 


d  v  _  du 

d  (  dc  + 


v  ( 


1  +  F2  ) 
Z_ 


du 

dc 


2 

F 

y 
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d0 

d< 


B 


u  £?  -  v  lH 
dt  •  d( 


u2  +  v2 


R- 


B 


ij 

d  € 


Then 


A*  =  1  sin  2  n' 

fdpi 

d0\ 

\  2y  ; 

L 

idcl 

■f  1 

3 

a4 

,  f  sin  pi  si|i  o'] 

n 

y  -  ; 

y  cos 


sin  2  U 

2  r 


(p; 


'  PB}  +  0A  6B 


JA 


+  xA 


tan  (0  +  u)  -  tan  e 
A  B 


A  € 


y  R  -  y  A+  x  A  tan  (  Q  +  fO  -  XB  *an  C  R  (  1  +  sin  2  c  R  ) 


x 


B 


B 


tan  {  0  +  jJL  )  -  tan  e 


B 


A  € 


1  +  sin  2  e 


y  =y  +(x-x  )tan(0+y) 

c  A  c  A  A 


€c  =  *B  +  *  « 


'm  =  o  W 

C  ^  CD  CO 


yc2  -  c  y^2  |  to  be  stored  in  table.  (  see  page  £,.7  ) 


The  above  equations  for  x^,  y  ,  and  A  €  do  not  apply  if 
(  0  +  H  )  =  90°  +  1  °.  See  page  42  . 

A 


Now,  iterate  on  T  to  calculate  the  remaining  properties  behind  the 


shock  at  point  C. 
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Iteration  on  T  : 


Initially,  set  T  =  1^ 


Then,  F  =  cot  e 
Y  c 


0! 


"2  » 


y  {  r  +  l  )  M" 

OO  C  00 


0  =  r  +  i 

c 


r  -i 

c 

y  -1 

CD 


r 

c 


OO 


M 


OO 


B  = 


.  1 


r  +  i 


a 


-2  r  +  i  +0\ 

'  c 


C  I 


T  +  1 

c 


EM 

2 

2'l 

1  +  F 

y 

+  Of 

1  +  F 

y  1 

-  a 


1  +  F2  I  +  B 


oo 


CD 

1  -(1  -  — ) 


2J  P  2  -1 

+  P  W  (  1  *  -2“  )  (  1  +  F  ) 

00  OO  00  a:  V 


h1 

c 


r  P 

c  c 

\1  1 

p 


r  -n 

h 


y  c  ; 


4. 506*10 

.  1 

p  =  2.302585 


•  l  » 


4n 


P1  =  i 

c 


n 


2116.4 


0. 002498 


Test  for  either  dissociated  gas  analysis,  page  38  ,  or  vibrational  excitation  gas 
analysis  page  39  . 
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Dissociated  Gas  Analysis  -  Detached  Shock  Point 


4  5 

I  E 

i=0  j=0 


Aij  (  P'  )'  (  h'  )J 


Ap' 


4  5 


A  T  =  -2.30259 


I 

m=0 

i 

n=0 

5 

I 

m  . 

n=0 

( 

ipj 

\ 

m  ,  n 


(m-1)  n 


(  r  +  i  )  ap ;  / 

C  C  \ 


+  — (  1  +  F2  ) 


r  +  i  r 

c  c 


B 


1 


+  (  1  +  F  ) 

( r  +  i  )2  y 


( i-n 


a 


1L 


(r+  i)  r  (r-i)  yu  (r  -  d 

C  C  c  CO  00  c 


a  1 


{  1  +  F  ) 

y 


-1 


r  =  r  +  at 

c  c 


Iterate  on  T  until 
c 


AT 


<  10 


Continued  on  page  40 
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VIBRATIONAL  EXCITATION  GAS  ANALYSIS  -  DETACHED  SHOCK  POINT 


y  - 1 

'  oo 


0'/T 


Iterate  for  Tc  ,  with  Tc  initially  = 


If  07T  >  70,  then  T  =  — ' 

J  c  c  C 


*P'c  = 


in  I — I — I  where  T  =  _ P  c_ 


2.30259 


AT  =  -  2.  30259 


( rc+  i)  aP'c 


r  +  i  r 

c  c 


(  1  +  F2  ) 


+  —  - - - 5  +  (l+^)LSl —  _U^£L+_2J +  1  U2  (  14-tt2) 

B  (  r  +  i)  r+  i  r  r  - 1  y  m2  (r  - 1)  r 

Lc  'C  C  C  GO  CD  C  C 


T  =  r  +  AT 

c 


Iterate  on  T  until  ATj  <  10 
Continued  on  page  40 
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The  final  values  of  P '  ,  p  ,  j)’  ,  p  ,  h'  ,  h  are  now  calculated,  using 

c  c  c  c  c  c 

either  dissociated  gas  analysis  or  vibrational  excitation  gas  analysis,  depending 

c 

upon  the  physics  at  the  particular  shock  point.  Note  that  ( — f  has  been 

R  c 

calculated  for  the  dissociated  gas  analysis  only;  and  has  been  stored  in  table. 
Now  we  may  compute  0^  : 


u 


W 


oo 


■  (  1 


00  N 


(  1  +  I^) 


W  F  (  1  - 
oo  y  ' 


00\ 


v 

c 


(  1  +  F  ) 
y 


0  =  tan 

c 


-1 


u 


It  is  now  necessary  to  calculate  the  parameter, 


de 


/  P  \ 
1  oo 

l  P  C  I 


This 


parameter  must  be  calculated  at  point  C,  at  this  time,  so  that  it  may  be  used 

during  the  course  of  calculating  the  next  point  on  the  detached  shock.  This 

parameter  is  required  to  obtain  (see  page  35  ). 

’  d  c 

This  is  done  by  differencing  this  point  C,  with  a  shock  point  with  an 
inclination,  -  1°  ,  point  2.  The  iteration  on  T  is  repeated  for  point  2  to  obtain 


p  .  Then 
2 


d 

1  - 

dc 

\PC 

(  1  - 


00 


)  -  ( 1  - 


oo 


A  € 


A  *2  =  '  1 


Then  compute  ,  y ^  ,  and  as  for  an  interior  point  (  gas 

dissociation  ),  page  45  ,  or  (vibration  excitation)  continue  on  page  41  . 
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VIBRATIONAL  EXCITATION  GAS  ANALYSIS  FOR  DETAC1&D  SHOCK  POINT 


where  Ec 


e0/4c 


This  value  of  entropy  is  to  be  stored  in  the  table. 


>  70 

I 


set  Ec=  1  and 


0' 

!f  >  70  ,  set  Ej,,  =  1 


In  either  case, 


In 


°°  -l' 
c  -1. 


il 

T»  "Tc 


W 

C 


acMc 


where  ac  is  found  as  for  an  interior  point 


.  -1 

pc  =  sin 
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SHOCK  POINT  CALCULATIONS  WHEN  (  6  +  u)  = 


x 


c 


X 


A 


Ac  t 


sin  2|a  \ 


pi  _pt 

A  B, 


(  V  °B>  '(yc  -  y  A> 


sin  (A  sin  0 

*  I A 


sin  2 

2y 


dP1 
d« 

A  ’  'B 


d£ 

de 


y 


c 


B‘ 


tan  e 


1  + 


A  e 

- - 

sin  2  e  g 


Note  that  in  the  expression  for  A  e,  the  last  term  in  the  numerator  contains 
yc.  As  a  first  approximation  compute  Ac  ,  omitting  this  term,  then  compute 

yc  ,  then  re-compute  Ae:,  using  this  value  of  yc. 


43 


Calculation  of  Interior  Point 


The  following  properties  are 

known  at  points  A  and  B: 

S 

x,  y  ,  P',  0,  R  ,  h'.p'.y,  (j,m  ,w. 
The  object  is  to  generate  these 
properties  at  point  C. 


y  B  ■  y  A"  XB  tan  ^  )b  +  XA  tan  (  ®  +  JLt  ^A 

x  = - 

C  tan  (  0  +  (X  )  -  tan  {  0  -  JX  ) 

’  A  B 


y  C  =  7  a  +  (  Xc.  “  XA}  tan  {  8+  ^  }A 


= 

1 

1  sin  2  (x  pi 

-4-  - 

sin  2  /x 

I  sin  2  fx 

I  +  1 

sin  2  /j  \ 

l  2  y 

1  ■  I 

A  ' 

2  y 

I  2y 

A 

1  2  r  B 

■+  eA  -  0 - 

A  B 


sin  JX  sin  0 


L  y  cos  (  0+jix) 


( xc -XA>  ~ 


sin  n  sin  0 


L  y  cos  (  0-fx  )  j 


(  Xcr  XB) 


1  sin  2  fjL 

sin  fi  sin  0 

\  2  y 

J 

«Pc  -  PA  >  - 
A. 

l  y  cos  (0+ji)  ^ 

The  above  equations  do  not  apply  if  (  0  +  |X  )  =90  4  1° 

.A 

For  the  equations  governing  this  condition,  see  page  48. 
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m. 


=  X/2  l™c.  +-  me  ) 


where 


mCl  =  mA  +  PAWA  sin>AA-  1/2  (YA+  Ye) 


( x c “  x  A) *  +(Yc*Ya)‘ 


/2 


mc2  =  mB  “  PB%sin:>B-  ^2  ^B+yc> 


(xc-xB)  +(yc-YB) 


l/2 


Enter  the  table  of  mass  flow  vs.  entropy  with  the  above  value  bf  and 

S 

extract  the  corresponding  value  of  (  R  ) 

k.  c 

Test  for  either  dissociated  gas  analysis,  page 45,  or  non-dissociated 
vibration  excitation  gas  analysis,  page  46. 


% 


Dissociated  G3.S  Analysis 


Let  h'*c  =  \  <  h’A  +  h‘B  >  and  P* 

4  5 


then  Ah' 


1 

2 


_S_ 

R 


Y  Y  Aij  (  P'  )  (  K  )y 

i=0  3=0 _ 


r  r-^.  ^  1 

Y  Y  j  *  Aij  (p«)  (h^ 

i=0  j=0 


1  w  S^1' 


Ap'  = 


i'll 


< “a +  V 


4  5_  m  n 

Bmn  (  p'  )  (  h*.  )  (  1  +  . _  )‘ 


rfcpO  n=0 


h'. 


(1  £ 


(m-1)  n-, 

m  •  Bmn  (  p'  )  (  h '  ) 


and  h'  =  h'^  +  Ah’  y  p^  =  p'^  +  Ap' 

c  c  c  c 


h  =  4.  506  *  I06  h' 
c  c 


p  =2116.4  e 


P  1 


2. 302585  p' 


p  =  0. 002498  e 
c 


2  /,  .  c 

a  =  (h  -  e  ) 


c  c  er 


where  e  =  (h 
c  c 


W 


Wc  =  a/K1  -  2hc  ,  Mc  =  — 


jU  =  sin  (  __L)  ,  y 

M  c 


ac  Pc 


Continue  on  page  47. 


Non-Dissociated  Vibrational  Excitation  Gas  Analysis 
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Let  p*  =  1/2(PA+PB).  p  =1/2  (pA+pB),T 


P 


R  p 


Then  Tc  =  T  +  AT 


where 


A  T=T 


* 


s\ 

Is' 

1  ^  |Too\  J 

Pc  ’ 

+  r  fn  -TT  +lr\ 

+  fn 

lRl" 

lR, 

•oo  00  \  T"'  }  ‘ 

Poo, 

O-’/T  1 

e  - 1 

e'/T 

e  '  A  oo  _  p 


+  0 


oo 


oo 


E 

7 


Foo  + 


-®4-|  E*  (E*  -1) 
T* 


>  0 '/  t  * 
;e'/T*  7i 


where  E*  =  e"  > 1  and  1^  =  'oo 

y  ~  1 

oo 


If  >7  0.  E*  =  1  and  fn 


0  /T: 


-1 


O'  it 

e  /Aa>  -l 


0'  _  0' 
T*”  ^oo 


Then  pc  =  ££_  and  p  = 


Pc  J  _  1  ln  I  PC 


RT, 


2.302585  0.002498 


hc  =  CPoo 


Tc  + 


roo(e9'  /Tc  -1) 


and  h'c  = 


4.  506  -  106 


7c 


=  1  + 


T'co  -1 


1+  <>00  _1) 


,0-  2  e  ^ 

T  I  Q'/Tc  1'2 


M, 


iL 

Tc 

2  _ 


M: 


sin 


7cTc 

-1 


=  Cp  .  Tc  and  yc  = 

y  co 

roo 

j  i 

1  .  _  \ 

Too  <Tt  -  Tc  )  +  »’ 

e'/Tt 

e'/TC  •  , 

m 

e  -1 

e  -1  /  J 

l Mc, 


.  Wc  =  J  Kj  -  2hc 


I 


An  improved  value  of  the  mass  flow  may  now  be  calculated. 


rn  = 
c 


where 


= 


m 

c. 

2 


V2  (  m  +  m  ) 


ci  cz 


r 


m  A  +  1/2  I  Pa  WA  Sin  M  A+  ^  Sin  M  ■ 


(  yA  +  yc} 


m 


B 


'  'll  PBWBSin^B+  PcWc  3intic) 


(  yn  +  y) 

±3  C 


This  completes  the  analysis  of  an  interior  point 
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o  o 

Interior  Point  Calculations  when  {  0  +  4)  =  90  +  1 


x  =  x 
c  A 


y 


c 


yB  +  ,\  'V  tan 


9-<*>b 


i 

1  sin  2 

JL  p,l 

sin  2  4  J 

P‘  =  1 

[  sin  2  (x\ 

I  sin  2  fi 

1  2  y 

j 

!  + 

A 

Jr' 

-  2y  jj 

c 

1  2  y  i 

+ 

'a 

\  Zy  1 

B 

1  sin  u  sin  0^ 

1  sin  4  sin  0 

-(yc- 

^a) 

l  ; 

y 

1 

1  -  ( 
A 

x  -  : 
c 

K  ) 

b' 

'  y  cos  (0-  4) 

/  sin  2  a) 

sin  4  sin  0 

c=0A- 

1  2  y 

(Pc-PA}  -  <yc  'V 

y 
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POINT  ON  BODY 


The  general  equation  of  second  degree 

.  2,2 

is:  ax  +  by  +  cx  +  dy  +  e  =  0 

This  equation  may  be  used  to  prescribe 
the  sequential  segments  which  make  up  the 
body  profile. 


■3 


2D _ 

-bp/bHad 


for  0  ; 


x 


A 


0 


Yc  =  Yb  +  (  xc *xg)  tan  (.  6  -n)B 


where  A  =  a  +  b  tan  (  0  -  p.) 


B 


B  =  tan 


(0 


yB"  xBtan  *  0 


+  d  c 


D  =  b 


yB-  xB  tan  (  0  -p.)  J 

2 

+  d 

Y  B~XB  tan  (  G  "^)  B 

>■  1 

• 

+  e 


0C  =  tan 


Pc  =  PB  + 


2  a  x  r  +  c 
2  b  y  +  d 


sinp,  sin  0 

Oc-Ob  ~ 

L.y  cos  (  0-p.) . 

sin  2  [i 


J 


B 


S 

R  A 


Test  for  either  dissociated  gas  analysis,  page 45,  or  vibrational  excitation  gas 
analysis,  page  4fi 


50 


EXPANSION  CORNER 

When  an  expansion  corner  is  detected,  the  program  adjusts  the  mesh  of 
characteristic  lines,  so  that  a  second-family  line  intersects  the  body  very' close 
to  the  corner.  The  exterior  angle,  a,  is  divided  into  15  components,  i.  e.  , 

A  Q  a 

~T5*  *  The  procedure  then  follows  a  step-wise  integration  around  the 
corner  to  form  the  vertex  of  the  expansion  fan  . 


Given  that  point  B  ife  the  ’'previous"  point  on  the  corner,  the  properties  of 
point  C,  the  " next"  point, are  found  as  follows: 


xc  =  x  B  > 


/  S  1 

f  s 

Yc  =  yB>  (r) 

C  I 

lR 

B 


©c  =  +  A  0  (Note  that  a  and  A  0  are  negative  fpr  an  expansion  corner.) 


Sr±_ i) 


Z  sin  Zp  cos  p 


-  sin  p 


Test  for  either  dissociated  gas  analysis  page  45  or  vibrational  excitation  gas 
analysis,  page46,  to  complete  the  calculation  of  the  remaining  properties  at 
point  C. 
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SECONDARY  SHOCK  CALCULATION 


When  the  program  detects  a  re-entrant  corner,  the  mesh  is  adjusted 
so  that  a  second  family  line  terminates  very: '.close  to: tH e  corner.  Refer  to 
Figure  Ilia,  page  95  .  Points  A  and  D  refer  to  the  upstream  and  the 

downstream  properties  at  the  corner,  respectively.  Point  A  is  calculated 
in  the  usual  manner,  utilizing  the  properties  of  points  F  and  J. 

To  cross  the  shock  at  the  corner,  it  is  necessary  to  estimate  the 
angle  of  inclination  of  the  shock  ,  Defining  F  =  cot*(  -  0^  )  ,  we 

solve  the  cubic, 

Fj>  +  tan  6  {  1  1-  M2  )  F2  +  ( 1-M  2  )  F  +  tan  6  (  1  +.yd,.T_lM 

v  2  v  n  n  2 


2 

U 


) 


where  the  subscripts  u,  d, pertain  to  upstream  and  downstream  properties  at' 

a  shock  point,  respectively.  Thus,  at  the  corner,  u  pertains  to  point  A, 

d  to  point  D.  This  cubic  implies  that  yd  =  yu  »  thus  it  yields  an  approximate 

value  foj:  F  .  Theflow'.deviation,  &  =  ct=  0„  -  0,  ■ 

»/  DA 

It  is  now  possible  to  cross  the  shock  to  obtain  the  downstream  properties 
at  the  corner.  This  analysis  is  similar  to  that  used  to  cross  the  detached 
shock  (page  37  )  except  that  upstream  conditions  must  be  substituted  for 

free  stream  properties.  The  iteration  is  performed  on  rather  than  on  T. 

At  the  corner  (only)  an  additional  iteration  is  performed  on  F  ,  since,  as 

V 

indicated  above,  the  cubic  yields  only  an  approximate  value. 
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Iteration  on 


Initially,  set 


a= 


Vd 


yu  (  yd  + 1)  mu2 


ft- 


B  =  ^ 


1  - 


/ 

1  1 

2 

a 

rd+ 1 1 

-  2 

?d  + 

2L\  =  - 

1 

-  /V 

yd  -  1  Pk  -  kk 

V+1  \V  yd-*/yuMu 

\ 


+  3 


l  +  F2'  |  +  a2 

V 


Wf*'  2 

tj 

\  / 


Pd  ?d+ 


r-' 

t? 


u 


Pd  1  -  (  1  ) 


pd 


/Xi  \ 


V1 


p-j 


P  =  p  +  p  w 

d  u  u  u 


d  =  4.506. 106 


1 


d  2.302585 


in 


Pd 


0. 002498 


d  ~in 


P 


d 


12116. 


i/, 


t 


1  - 


P;d 


1  +  F 


Teit  for  either  dissociated  gas  analysis  or  vibrational  excitation  gas 

analysis. 

Dissociated  Gas  Analysis  ^ 
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m 


n 


6  p' 
d 


I 


z  i 

m=0  n=0 


Bmn 


4 

5 

— 

m-1 

n 

I 

m=0 

\  m  *  Bmn  Ip' 

L -*  d 

n=0  ' 

h1 

d 

y 

>  (2a) 


Continue  with  equ.  (3) 

Vibrational  Excitation  Gas  Analysis 


T  = 
d 


I  y  -1 
\  'co 


0' 


CO 


Pd 

P  R 
d 


^oo 


Apd' 


(  0'/T) 
e  -  1: 


1 


Iterate  for  T 


2. 302585 


In 


After  obtaining  Ap'  using  either  (2a)  or (2b), 
d 


Ayd  =  -Kc 


P  l  a 

Pd  d  d  |  jyd  + l)  yd 


>  (2b) 


(1+pf)  +  ± 

7?  B 


_ I -  +  (  1  +  FZ 

L(yd+  n*  v 


-1  (3) 


a 


( 1  -yd) 

- +  -  + 


a' 


t(y+ l)  y, 

d  d 


(  1  +  F2  ) 

y  - 1  y  M2  (y  - 1)1  y  17 

*  u  u  d 


,2  *]  V 


Then  y  =  y  +  Ay. 

d  d  d 


K  =  2.30259 

c 


The  iteration  on  v  converges  when 
d 


Ay, 


<  10 


-6 


For  the  corner  only,  equ.  (4)  is  executed;  for  all  other  secondary 
shock  points,  the  analysis  continues  with  equ.  (5). 
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AF  =  J  tan  6  [  {  1  +  F2  ')  -  -Cl  (  1  +  F2  ' )  +  b]  1-  — 2- 

V  |  L  V  yd  +  1  n  V  1 


+  a(  1  +  F2  )  F  -BF  >  * 
V  V  »» 


<  rVi  “ ( 1  +  ^ )  +  B  -2F>! ta” 8  +  21 V  *an8:fe( “( 1  +  pn>  *r«i]j 


F  =  F  +  AF 

V  V  V 

The  iteration  continues  by  returning  to  equ.  (1),  and  converges  when  A  F  J  <  10 


u=  Wu  1 


( i  -  ^ ;) 


(  i  +  f  2: 


(  i  -_h  ) 

v  =  WuFv  Pd 

(  1  +  F*  ■ ) 


- 1  v 


0,  =  tan  (  u  ) 


Test  for  either  dissociation  gas  analysis,  equ.  (6a)  ,  or  vibrational  excitation 
gas  analysis,  equ.  (6b). 


Wd  =  VKi-2hd 


M  d  =  — - 
d  a 


•  -1  L-\ 

^d  =  sin  (  Md 


where  a2  =  (  h  d  -  e^)  _ 


ed  =  hd 


y*.  a  Pd 
— - 

pd 


[ 
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L  In  either  case,  y ^  is  written  as  output  and  used  in  the  calculation  of 

[ 

[ 

i 

La 

[ 


properties  of  points  downstream  of  the  secondary  shock 
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Calculation  of  Point  on  Secondary  Shock  above  Corner 

After  crossing  the  secondary  shock  at  the  re-entrant  corner,  a  new  second 
family  line  is  generated,  originating  at  the  detached  shock  and  terminating  at 
point  C  (  see  Fig.  Ilia)  ,  Note  that  A  -  C  is  a  first  family  line  which  is  a 
part  of  the  mesh  upstream  of  the  secondary  shock,  i.  e.  the  properties  at 
points  A  and  C  are  those  which  would  exist  if  the  secondary  shock  were  not 
present.  Then, 


y  -  x  tan  e  -  y  +  x  tan  (  0  -  u  ) 

A  A  o  B  B  '  p  B 


tan  (  0  -  n  )  -  tan  e 

B  ° 


Y 


B' 


y  +  (  x  -  x 
B  B'  B 


tan  (  0  -  n  ) 
B 


The  upstream  properties  at  B'  are  found  by  interpolating  between  points 
B  and  C. 

Using  the  downstream  properties  at  the  corner,  0  ,  and  a  first 

family  is  generated  connecting  B'  with  the  body,  at  point  G. 


2D 

x  =  - — rzrzzzur' 

G  -  B  +  Vb2  -  4  AD 


for  A  ^  0;  x  = 
G 


for  A  =  0 


yG=yB'+(XG‘  Vtan  {9  +  MD 
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where 


A  =  a  +  b  tan  2  (  0  +  ju 
B  =  tan  (  8  +  n  f  2  b 


D  =  b 


i'B i  -  XB>  tan  (  9  +  li 


yB,  -  xB.  tan  ( 0  +  m  )D 


+  d 


+  d  )  +  c 


yB,  -XB-  tan  (0  +  fi  )d'+6 


2  a 

x  +  c 

“ 

G 

2b 

y„  +  d 

SG  •  tan 


Note  that  the  coefficients,  a,  b,  c,  d,  e,  refer  to  the  body  profile 

immediately  downstream  of  the  corner,  A-G-E. 

The  remaining  properties  at  point  G  are  assumed  equal  to  those 

properties  at  point  D. 

Point  ^  is  located  by  the  intersection  of  a  line  drawn  from  the  corner 
having  an  angular  inclination  of  (  c  -f  -y-  )  ,  and  the  second  family  line  joining 
points  B  and  C.  The  value  of  A  is  chosen  to  be  -1°»  Then, 


x~ 

B 


YA  ~  XA  tan  (  €o  }  ~  yB+  ^tan  {  6~  fi) 
tan  {  0  -  jU  )g  -  tan  (  co  +-y  ) 


y0  +  (  x«  -  x  )  tan  (  0  -  M  )n 


B 


B  B 

l 


B  and  C , 


The  upstream  properties  at  B^  are  found  by  interpolating  between  point 

Point  Bx  is  located  by  the  intersection  of  the  first  family  line,  B'  G,  and 

A  ~ 

the  line  having  the  inclination,  (  c  +  -g-  )  ,  AB^  . 


% 


yA  -  XA  *an  <  6o+  -f  >  ~  yG  +  XG  tan  (  6  +  M  }G 

tan  {  Q  +  n  )  -  tan  (  c  +X  ) 

G 


yBj  =  VG  +  (  XB,'  *G  1  *“  '  9+  "  ’g 


58 


The  upstream  properties  at  B[  are  found  by  interpolating  between  Bx  and  A. 

If  the  secondary  shock  were  straight  from  the  corner  to  its  intersection  with 
the  first  family  line  through  G,  then  B"  would  be  a  point  on  the  secondary  shock. 

On  the  other  hand,  if  the  secondary  shock  changed  in  curvature  by  the  amount  A  . 
between  the  corner  and  its  intersection  with  the  first  family  line  through  G,  then 
Bj  would  be  a  point  on  the  secondary  shock.  The  actual  intersection  of  the 
secondary  shock  and  the  aforementioned  first  family  line,point  B,  lies  somewhere 
in  the  neighborhood  of  these  two  points,  B1  and  Bj  , 

To  locate  point  B  ,  and  to  calculate  the  downstream  properties  at  this  point, 
the  following  procedure  is  followed. 

1.  Assume  that  B'  lies  on  the  shock,  and  "cross  the  shock"  to  get  downstream 
properties  at  B' 

2.  Assume  that  Bj  lies  on  the  shock,  and  "cross  the  shock"  to  get  downstream 
properties  at  Bt 

3.  With  this  information,  the  location  and  properties  of  the  actual  point  on 
the  shock,  B,  may  be  calculated. 

To  get  the  downstream  properties  at  point  B1,  set  F^  =  cot  (  -  0gf)  ,  and 

iterate  on  y^  (  equ.  1-3),  omitting  equation  (  4)  since  F^  is  prescribed.  The. 
remaining  downstream  properties  are  found  by  equations  (  5-6)  . 

Similarly,  to  get  the  downstream  properties  at  point  Bj  ,  set  F  =  cot'(  c  +  A  -0  ) 

V  °  Bj  ' 

and  proceed  as  above. 


Now, 


B, 


-  P 


dP' 
d  c 

dh' 
d  ( 


B' 

d 


de 

d  € 


^d  ’  9B'd 


(7) 
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The  actual  change  in  curvature  between  points  D'  and  B  is 


Ac  = 


sin  2  ji\ 

iTy~iG 


Pc-%n-w 


sin  sin  9 
T  cos(9fji)jQ 


x„  +  Y  G~  xctan  (  9+  M)g-  y  D+3!Dtan<b 
G 


tan  (  0  +  H)r-  tan  e 


sin  2  u 

f 

d  P 

2  y 

G  d  C 

d  0 
d  € 


-  1 


A  e 


X- 

B 


y  t)  ‘  XD^  tan  €0  ^  (  1  +  sin  2  )  -  y.  +  xG  tan  (  0  +  JX  )  Q 


A  € 


tan  (  0  +  n  )  -  tan  c  (  1  + - - — - - 

G  ^  sin  2  Co 

yB  =  yG+  (XB  -XG  )tan(  0+^  )Q 


) 


dP 

PB  =PB  +  'dt  I  *  ‘ 


I  I 

h  -  :  h  + 

B  B 

d 


*S  =  €o  +  A€ 


dh 

\ 


A  e 


9B 


=  0  + 
B'd 


d_0  \ 
d 


A  € 


(  7)  cont. 


(8) 


Test  for  dissociation  gas  analysis,  or  for  vibrational  excitation  gas 
analysis. 

For  dissociation  gas  analysis,  use  equs.  (  2a  ),  and  (  6a  )  to  generate 
the  remaining  properties  at  point  B.  For  vibrational  excitation  gas  analysis,  use 
equs.  (2b),  and  (  6b  )  to  generate  the  remaining  properties  of  point  B. 

Note  that  these  are  downstream  properties  at  B, 

After  computing  the  downstream  properties  of  B  ,  the  body  point,  E,  down¬ 
stream  of  the  secondary  shock  is  calculated  in  the  usual  manner,  page  49  .  In 
t  his  case,  point  B  acts  as  point  B;  and  point  D  as  point  A(page  49  . ).  The  line 
B  E  is  a  second  family  line. 

The  program  then  returns  to  the  detached  shock,  to  generate  a  new  (upstream) 


second  family  line. 
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CALCULATION  OF  SUBSEQUENT  POINTS  ON  THE  SECONDARY  SHOCK 


Refering  to  Fig.  Illb,  the  point  D  on  the  secondary  shock  has  been  computed, 
as  well  as  point  E  on  the  body.  The  second  family  line  terminating  at  point  C  has 
been  generated. 

To  determine  point  B,  a  procedure  similar  to  that  used  to  calculate  point 
D  is  followed.  With  cq  =  point  B1  is  calculated  in  a  manner  identical  to 

that  described  lot  the  shock  point  D,  above  the  corner. 

Point  G  is  the  intersection  of  the  first  family  line  from  B'  (using  0  ^  and 
/Lt  £)  ),  and  the  downstream  second  family  line,  DE. 


Y  D  "  y  B  "  x  D  tan  (  9  “  M  )  D  +  XB'  tan  (  8  +  ^  )d 

tan  (  0  +  ii  )  -  tan  (  0  -  ) 

D  D 


Y  , +  (x  -  x  )  ton  (  0  +  H  ) 

B  G  Bf  D 


The  remaining  properties  at  G  are  interpolated  between  points  D  and  E, 

The  points  B^  B^,  and  B  are  calculated  in  exactly  the  same  manner  as 
described  for  the  shock  point  above  the  corner,  pages  56  to  59  . 

The  second  family  line  downstream  of  the  secondary  shock,  is  now  generated. 
The  interior  point,  F,  is  calculated  using  points  B  and  E,  in  the  usual  manner,  where 
B  acts  as  point  B,  and  E  acts  as  point  A  (page  43).  The  body  point,  H,  is  found 
using  point  F  and  E. 

In  this  manner,  the  secondary  shock  is  constructed  along  with  the  rest  of  the 
flow  field.  Note  that  the  second  family  lines  downstream  of  the  secondary  shock, 
are  slightly  offset  with  respect' to  the  upstream  second  family  lines. 
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When  there  is  only  1  interior  point  on  the  second-family  line  downstream  of  the 

first  secondary  shock  then  we  must  test  whether  point  Q  lies  below  point  P  or 

Test  y  ^  y 

Q  P 


where 


yA  -  xA  tan  CA  -  YB  +  tan  C] 


tan  €  -  tan  e 

B  A 


;  yQ  =  yB  +  (  *Q-  XB> tan  € 


B 


I 


yA  ~  XAtan  (A~  yE+  Van  (  9~  ^  }E 
tan  (  6  -  H  )  E  -  tan  *A 


•  yP=  yE+  (  xp_  xe)  tan  (0'M)i 


If  y  •>  y_  ,  no  intersection  takes  place  and  the  2  s.  s.  must  be  extended  at 
'Q  P  c 

least  one  more  point. 


I 


If  y  <  y_  ,  then  the  shock  intersect  as  shown. 
Q  P 
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The  conditions  upstream  at  point  Q  is  assumed  to  be  those  upstream  at  point  A 


(stored  in  PNTF,  which  is  B1)  :  W_  =  W_ 

Qu  F 


It  is  necessary  to  get  an  approximate  value  of  the  inclination  of  the  resultant 
secondary  shock  of  Q: 


,  rQ  +1  \  /  r  q-1  ' 

F  3  +  tan  6  1  +  u  M 2  F2  +  (  1  —  M  2  )  F  +  tan  6  1  +  _J± - M2  , 

V  o  \  2  QJ  V  Qu  V  °  \  2  Quj 


=  0 


where  6  =0-0  initially 

o  B  F 


Then 

-1 

€  =  0  +  cot  F 

f  n 

Subscripts:  Q  ,  downstream  at  point  Q 

d 

Qu  ,  upstream  at  point  Q 
B  ,  downstream  at  point  B 
A  ,  downstream  at  point  A 
F  ,  upstream  at  point  A 
n  ,  for  n  n  iteration  on  6 


(1) 
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Now,  for  the  given  value  of  6  ,  we  may  cross  the  shock  in  the  usual  manner: 


'Q 


oc  = 


1 


yQ 


U 


(  y  +1)  M2 
~u  00 


p  - 


v<V‘  lY° 


u 


y^  +  1\y^  -l  y^  -  1  M 


Q 


u 


Q  Q 
Mu  u 


Pn  \ 


Q 

1  -  - ^ 


/ 


— - -  -  a  (  l  +  F“  )  + 

\+1 


1  -0(1+F*) 


LyQd  +  1 


-2(3  (  1  +  Fj)) 


Q 


u 


Q, 


1  -  1  - 


V 


Qd 


(2) 


PQd=  PQU+  W 


1  - 


(  1  +  F2  ) 

V 


-l 


Q 


yO  _1 
ud 


Now  test  Real  gas  (3a)  or  vibration  excitation  (3b) 

P, 


1 


2. 302585 
h 


QH 

in  -  _d-_  _  ,  P;  ±  in 


Q, 


hQ  ■ 


Q. 


2.498. 10-3 
I  S  \  _ 


4.506-  106 


R 


d 

4,  5 

£ 


2.  1 1 64 r  i 0 3 


*p: 


Q, 


AU  (PQd)1(hbd)J  )(3a) 


64 


T  = 


Q, 


co 


yoo  -i 


'00 


eVT 

e  - 1 


by  iteration  T.in 


"Poo 


4=1; 

ud 


y 

GO 


y  -1 

00  00 


'in 


;t  ^ 

T 

>  oq I 


+  in 


lp  00  ' 


+  9' 


1 


2. 302585 


\pqJ 


in 


P '  =  in 

Qd 


2. 1164-103 


0'./T 


0'/T 
e  oo 


+  in 


0’ /Too 


0'AT  . 

e  - 1 


(3b) 


Pod  1 

1 

—  1  _  in 

P°d  \ 

’  Qd  2.302585 

i2. 498- 10" 3 

h'  =  - 3 - 7- 

Qd  4.506-10° 


Then  after  using  either  (3a)  or  (3b)  , 


Ay^  =  -  2.302585  (  y„  +  1)  (  Qu /PQ_  > 


Q. 


'Q 


- I - +  _9L_  (  1  +  F2  ) 

+ 1  v 


i 

+  B 


where 


1  -7t 


'  Q  d 


M+  A 


rQd  ’V 


2  + 


VqjJ 


--21  |  i^J* 

■Vr. .  » 


B2  = 


-2 


a  \ 


(  y„  +  i  J2  lr_  +  1 

^d  / 


+  0  (  l  +  F2  )  +  a2(  l  +  F2  ) 

r?  V 


and  y  =  y  +  Ay 

Q,  Qa  Qj 

d  d  d 


ft  4) 


With  this  new  value  of  y0  ,  repeat.eqs  .  (2)  thru  (4)  until  Ayn  s  i0 

Uj  MJ 


With  the  final  iterated  value  of  y_  ,  calculate  the  correction  to  F  corresponding 

Qd  T< 


to  this  value  of  5: 
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tan  fi  l  (  1  +  Fp  - 


AF 


— - —  -a  (  i  +  r2)  +  b 
y  +  1  V 

rD 


tl 

- —  +&  (  1  +  F2  )  F  -BF 

y  +1  T?  T?  T) 


V  - i —  -  a  (1+F2  )  +B  -  2F  tanfi  +  2(F  +tan6) 

y  +  l  t)  r? 


-^(c^  1  +  F2  )-S-— )-«F 
B  1?  y  +1 

D 


and  F  =  F  +  AF 
V  V  7) 

With  this  new  value  of  F  repeat  eqg.j.  (2)  thru  (5),  iterating  on  both 

V. 

AF  Is  10'5 

V 


and  F  as  indicated  until 


(5) 


With  the  fiqal  iterated  values  of  y^  and  F^ ,  test  for  real  gas  or  vibration 
excitation  and  compute  4  ,  accordingly,  as  usual. 

Then  compute  a  new  valuer  for  6  : 


6„+l  =  «o 


sin  4  cos  H 


% 


Z)t  _  pi 

Qd  b; 


(6) 


n 


With  this  new  value  of  6  =  6  ,  ,  test  below: 

n  n+ 1 


6  -  6  + 
n  0 

> 

sin  n  cos  n) 

1  PI  .  Pl|  I 

I  r  , 

ki  1  B) ' 

n _ 

6n 

<  0.0025 


If  the  test  fails  (  >  0.0025)  return  to  (1),  etc. 

If  test  is  satisfied,  the  properties  of  the  resultant  shock  at  Q,  downstream, 
have  been  calculated.  Use  (1)  to  calculate  inclination  of  this  shock  of  Q 
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Intersection  of  Bow  Shock  with  Secondary  Shock 

The  analysis  is  identical  with  that  outlined  for  2  secondary  shocks  except 
that  eqv..  (1)  becomes 

€  =  cot  ^  F  =  cot  1  F 

Q  17  y 


and  the  subscripts,  Q  and  F  are  replaced  by  cd  which  indicates  free  stream 


u 


conditions,  and  0=0  =  o.  Thus  6=0-,  initially. 

F  o  B 
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CALCULATION  OF  ENTROPY  WITHIN  THE  FLOW  FIELD. 

As  indicated  in  the  preceding  analysis,  the  mass  flow  at  each  point  in  the 
mesh  of  characteristics  is  calculated.  The  entropy  is  explicitly  calculated  only 
at  points  along  the  detached  shock  and  along  all  secondary  shocks.  At  each  shock 
point,  the  calculated  value  of  entropy  and  the  corresponding  mass  flow  are  entered 
into  a  table.  Since  the  relationship  between  mass  flow  and  entropy  is  different 
behind  a  secondary  shock,  than  it  is  behind  the  detached  shock,  up  to  three  sets 
of  mass  flow-entropy  tables  must  be  constructed  -  one  for  points  downstream 
of  the  detached  shock  and  upstream  of  the  first  secondary  shock;  one  for  those 
points  between  the  first  and  second  secondary  shocks;  one  for  those  points  down¬ 
stream  of  the  second  secondary  shock. 

These  tables  are  logically  adjusted  and  modified  every  time  a  change  in  the 
character  of  the  flow  field  takes  place:  secondary  shocks  dying,  intersecting, 
originating  at  a  corner,  or  intersecting  the  detached  shock.  After  the  mass  flow  is 
calculated  at  each  interior  point,  the  program  correlates  the  position  of  this  point 
within  the  flow  field,  with  its  correct  entropy  table,  and  extracts  the  corresponding 
value  of  entropy  by  interpolation  within  this  table.  Since  the  mass  flow  is 
calculated  iteratively  at  the  interior  points,  these  values  are  of  high  precision, 
and  thus  the  values  of  entropy  extracted  from  the  tables  are  correspondingly  very 


accurate. 
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CONICAL  FLOW  ANALYSIS 

Given  the  free  stream  parameters  and  the  body  equation  for  the  conical 
nose'j  the  attached  conical  shock  angle  is  obtained  by  interation,  and  a  horizontal 
reference  line  consisting  of  20  points  is  generated  from  the  shock  to  the  body. 

From  the  body  equation  cx  +  dy  +  e  =0,  the  body  angle  is  obtained. 

TJ  =  tan  ’  *  (  -  /d  )  (1) 

c 

An  initial  guess  of  the  shock  angle  is  given  by 

n .  ...  ,  =  1.2  (TJ )  (2) 

initial  c 

and  r  is  assumed  =1.25 


Iteration  for  Shock  Angle 

With  r  and  F  =  cot  TJ  properties  behind  the  shock  are  obtained 
using  the  equations  given  on  pages  through  .  The  conical  flow  equations 

(7)  through  (15)  below  are  used  to  obtain  the  value  of  a  body  angle  (  TJ^  ).  for  the 
given  shock  angle  (  TJ.  ). 

For  the  first  iteration  on  the  shock  angle 


6. 

i 


=  T7. 

i+l  i 


+  6. 
i 


(3) 
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The  iteration  is  repeated  with  t?  =  t? 

i  i+ 1 


On  all  subsequent  iterations, 


s.  =  n  -n 

i  c 


A1  =  %  -  V  , 


(4) 


and 


EDELTA  =  {  6.)  (  fi.  j)  /  Aj 


7?.  =  T?.  +  DDELTA 

i+ 1  i 


(5) 


The  calculation  is  repeated  until 

IddeltaI  <:  .01° 


(6) 


For  the  final  shock  angle  TJs  =  TJ.  +  DDELTA  the  properties  behind 


the  shock  are  obtained. 


Conical  Flow  Equations 


Given  T?c>  the  shock  angle  t?  and  the  properties  behind  the  shock 


At?  =  (  rj,  -  7?  )  /  10 


(7) 


If  At?  S  1°,  then  At?  =  l.°  and  DETA  = 


At?/  100 


r  ,  ,  1 

—  =  {  u  cos  T?  +  v  sin  T?  )  - 

'  W 

00  co 


(8) 


70 


At 


=  (  -  u  sin  TJ  +  v  cos  r) 


(9) 


and  an  iterative  procedure  is  used  to  obtain 


starting  with  the  shock  angle  77 


from  the  following  equations 


V  \ 

/ 

_t  1 

- 

W 

" 

\  °°i 

—  cos  At)  +  —  -  —  sin  A  TJ 


W 


IVi-1 


co 


W_ 


co 


- 1 


(10) 


with 


w 

lw 

\°°i 

H 

sin&n 
Vi  1 


cos  A  77  + 
'j-1 


IR\ 

W 

'S-i 


1  vt  vr  \ 

w 

\°°i 

vT  cotTJ+w 
^  '00  **00  J 

1  _ 

VVWco 

z 

a/w«, 

z 

(11) 


(12) 


\ 


/hi 

- 

v  \ 

z 

V  \ 

2  -v 

1 

1  - 

r 

t 

W7" 

2 

W 

W 

\  00  / 

_ 

\  00  / 

\  °°i 

00 


(13) 


(14) 
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when 


0 


A  tj  =  A  tj  /  10  until 
|  A  T7 1  £  DETA 


V. 


is  then  interpolated  for 


(15) 
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Points  along  the  Reference  Line 


Once  the  final  shock  angle  (  T?  )  and  the  properties  behind  the  shock 

s 

are  obtained,  a  horizontal  reference  line  consisting  of  20  points  is  calculated, 
r  and  /R  are  maintained  constant  and  equal  to  the  shock  point  values.  The 
x-  coordinate  of  the  last  point  on  the  reference  line  (body  point)  is  4  of  the 


distance  to  ^  . 

Thus,  x(20)  =  —  (  Xj  -  ^  )  (16) 

y  (J)  =[(c)  x  (20)  +e]  /d  1  £  j  <;  20  (17) 

x  (1)  =  cot  T?g  y  ( 1)  +  xq  (18) 

Ax  =  4/  71  [x  (20)  -  x(l)]  (19) 

x  (2  ')  =  x  (1)  +  Ax  /  4  (20) 

x(  3)  =  x  (2)  +  Ax/2  (21) 

x(.J)  =  x(J-l)  +Ax  4SJ  s  19  (22) 

Values  of  17  (J)  corresponding  to  x  (J)  can  be  obtained  from 

T)  (J)  -  tan-1  fy(J)  /  (x  (20)  -  x  )  1  (23) 

L  O  J 


The  values  at  each  subsequent  point  J  are  obtained  in  the  following  manner. 

At?  =  [t?(J)  -  T?  (J-l)  ]  /  5  (24) 

If  At?  *  .01°  ,  AT?  =  At?  /2 

Equations  (10)  through  (14)  are  iterated  upon  from  r?  (J-l)  until  t?(J) 

is  reached  and  the  values  of  h/W*  ,  v.  /  W  .  v  /  W  ,  a2  /  W2  obtained,  at  each 

00  t  cxj  r  00  00 


point  (J) 
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u  = 


v  = 


0  = 


V 

V, 

r 

w 

oo 

cos  T)  - 

t 

sin  TJ 

1 

OO 

(25) 

V 

Vf 

r 

W~ 

sin  t)  + 

1 

w" 

cos  T? 

w 

00 

(26) 

l°° 

00 

-1 

tan 

('  V/u  ) 

(27) 

M 


(28) 

(29) 


Use  either  (30a)  for  gas  dissociation  or  (30b)  for  vibrational  excitation. 


Assume  p'  (J)  =  p'  (J  -  1)  initially  to  obtain 


/R 


I 


6  p' 


m=0  n=T 


Bmn  (fi.)  ( h')n 

1 


m  Bmn  (p.1  f  (h') 


m  =  0  n=t 


Iterate  using  p|+  ^  =  q}  +  Ap' 
then  obtain 


p  =  ph 


!>  (30a) 


r 
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INPUT  FORMATS 

/ 

In  this  section,  the  input  formats  for  each  of  the  three  program  decks  - 
Flow  Field,  Transonic-Subsonic,  and  Supersonic  -  will  be  described  in  detail. 
Refer  to  the  section  on  Nomenclature  for  allied  information. 

The  term  "card",  refers  to  the  standard  IBM  data  processing  card 
consisting  of  12  rows  and  80  columns.  Since  columns  73  -  80  are  not  read  by 
the  computer,  any  identifying  information  may  be  punched  in  these  columns.  The 
term,  "format",  refers  to  the  mode  of  input.  Symbolically,  these  modes  may 
be  defined  as  follows: 

I  integer  .  ±  XX  (no  decimal  point) 

F  fixed  point  ±  XX.  XXX  (decimal  point  required) 

E  floating  point  ±  X.XXXiYY  (YY  is  the  exponent  to  the 

(YY) 

base  10:  ±X, XXX-  10  ) 

For  the  E  and  F  modes,  the  decimal  point  may  ,  of  course,  be  shifted 
from  the  position  indicated  in  the  above  examples  and  the  maximum  number  of 
significant  figures  is  governed  by  the  field  width  assigned  for  each  "word"  of  data. 
The  plus  (+)  sign  may  be  omitted  in  all  cases,  except  for  the  sign  immediately 
preceeding  the  exponent  for  the  E  mode.  An  additional  format  is  the  Hollerith 
mode  which  consists  of  alpha-numerical  information,  and  for  our  purposes,  will 
be  utilized  exclusively  for  an  identification  input  card,  which  will  subsequently 
be  printed  as  a  title  at  the  head  of  the  output  listing. 
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FLOW  FIELD  PROGRAM 


Due  to  the  size  of  this  binary  deck,  an  option  has  been  incorporated 
into  the  coding  to  permit  the  user  to  either  read  in  the  input  data  cards  via  the 
IB  Monitor  input  tape  or  via  the  on-line  card  reader.  It  is  recommended  that 
the  latter  option  be  utilized,  since  this  will  allow  the  user  to  permanently  store 
the  object  deck  on  magnetic  tape,  thus  avoiding  frequent  handling  of  the  cards. 
Refer  to  the  section  on  Operating  Instructions  for  further  details. 


CARD 

COL. 

DATA 

FORMAT 

1 

1 

1 

I 

2-72 

Hollerith  information  -  title 

2 

1-10 

M 

CO 

F 

11-20 

T 

oo 

F 

21-30 

( 

F 

31-40 

*co 

E 

41-50 

Pco 

E 

3 

1-10 

X 

F 

11-20 

* 

y 

F 

21-30 

X 

c 

F 

31-40 

yc 

F 

41-50 

9c 

F 

51-60 

R 

c 

F 
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CARD 
4,  5,  .  .  .  ,  A 


A  +  1 


A+  2 ,  A+  3 


Col. 

DATA 

_ 

FORMAT 

2-10 

a 

F 

12-20 

b 

F 

22-30 

c 

F 

32-40 

d 

F 

42-50 

e 

F 

52-60 

Xt 

F 

62-70 

a 

F 

1 

Stagnation  Point 

Code  =  1 

I 

2-10 

X 

o 

F 

11 

Point  1  Code  = 

1  or  2 

I 

12-20 

X1 

F 

21 

Point  2  Code  = 

1  or  2 

I 

22-30 

X2 

F 

31 

Point  3  Code  = 

1  or  2 

I 

32-40 

X3 

F 

41 

Point  4  Code  = 

1  or  2 

I 

42-50 

X4 

F 

51 

Point  5  Code  = 

1  or  2 

I  ' 

52-60 

X5 

F 

61 

Point  6  Code  = 

1  or  2 

I 

62-70 

x6 

F 

Identical 

in  format  to  card  A  +  1 

(if  required) 


Pertains  to  point  s7  through  20,  in  sequence 
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CARD 

COL. 

DATA 

FORMAT 

B 

2-10 

Yo 

F 

12-20 

yi 

F 

22-30 

^2 

F 

32-40 

y3 

F 

42-50 

^4 

F 

52-60 

^5 

F 

62-70 

y6 

F 

B+l.B+2 

Identical  in  format  to  card  B. 

(if  required) 

Pertains  to  points  7  through  20,  in  sequence 

C 

1,2 

Number  of  intervals  of  mesh 

spacing,  A 

I 

3-10 

At 

i 

F 

Hi  12 

Number  of  intervals  of  mesh 

spacing,  AT2 

I 

13-20 

At2 

F 

21,  22 

Number  of  intervals  of  mesh 

spacing,  At3 

I 

23-30 

AT3 

F 

41,  42 

Number  of  intervals  of  mesh 

spacing,  Ay1 

I 

43-50 

Ayi 

F 

51,  52 

Number  of  intervals  of  mesh 

spacing,  Ay2 

I 

53-60 

*y2 

F 

61,62 

Number  of  intervals  of  mesh 

spacing,  Ay3 

I 

63-70 

Ay3 

F 
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Card  1  may  be  used  for  a  title  of  identification  since  it  is  written  as 
the  first  line  of  output.  The  free  stream  conditions  are  punched  on  card  2. 

Card  3  contains  the  coordinates,  (  x' ,  y  ),  of  the  assumed  sonic  point 
on  the  body,  i.  e.  where  the  slope  of  the  body  is  unity.  The  remaining  data 
consists  of  coordinates  of  a  downstream  point,  (  xc>  y^  ),  where  the  slope  of  the 
body  is  tan  0^  ,  and  where  the  radius  of  curvature  of  the  body  is  R_  .  In 
general,  20°  <  0^  <  25°  ;  the  input  value  of  0c  is  in  radians. 

Cards  4  through  A  contain  the  coefficients  of  the  equations  describing 

2  2 

sequential  sections  of  the  body  profile:  ax  +by  +  cx+dy+e  =  0.  The 
domain  of  validity  for  each  section  is  defined  by  its  terminal  value,  .  The 
exterior  angle  between  this  section  and  the  following  one  is  a  ;  a  positive  non¬ 
zero  value  of  a  defines  a  reentrant  corner  at  x  =  x^_  .  If  a  is  negative,  an 
expansion  corner  is  defined  and,  of  course,  a.  -  0  indicates  that  the  body  is 
continuous  at  x  =  x^_  .  A  maximum  of  12  body  profiles  may  be  prescribed. 

The  following  cards  prescribed  the  coordinates  of  up  to  21  points  which 
define  the  geometry  of  the  nose  region  of  the  body,  and  also  indicate  which  of 
these  points  are  to  be  satisfied  by  the  iterative  procedure  employed  by  the 
Transonic-Subsonic  program  to  achieve  an  accurate  configuration  of  the  bow 
shock.  A  maximum  of  10  of  these  "iteration  points",  including  the  stagnation 
point,  may  be  prescribed. 

The  x- coordinates  of  these  points,  in  sequential  order,  from  the  axis 
of  symmetry,  downstream,  are  punched  on  card  A  (and  cards  A+ 1  and  A+2,  if 
necessary),  along  with  codes  which  indicate  whether  a  specific  point  is  to  be 
satisfied  by  iteration,  A  code  =  1  prescribes  an  "iteration  point";  a  code  of  2 
indicates  that  this  point  will  be  used  for  the  purpose  of  describing  the  nose  geometry 
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only.  The  corresponding  y- coordinates  are  on  card  (s)  B  (  and  B+l  and  B+2  ) 
Note  that  the  stagnation  point  must  always  be  prescribed  as  an  "iteration  point". 

To  assure  an  accurate  solution,  it  is  advisable  to  prescribe  at  least 
two  points  adjoining  the  assumed  Sonic  point  on  the  body,  as  "iteration"points. 
For  a  spherical  nose,  these  points  should  be  located  where  the  slope  of  the 
body  is  approximately  40°  and  50°  respectively.  Should  the  body  exhibit  a  more 
rapid  change  of  curvature  than  does  a  sphere,  in  the  neighborhood  of  the 
sonic  point,  then  additional  iteration  points  should  be  clustered  in  this  region. 

At  least  5  points,  including  the  sonic  point,  should  be  prescribed  as  "iteration" 
body  points  in  the  region  from  the  axis  of  symmetry  to  the  assumed  sonic  point. 

The  final  input  card,  C  ,  prescribes  the  y -y  mesh  used  in  the 
subsonic  region.  As  indicated,  up  to  3  different  mesh  intervals  in  each  direction 
may  be  prescribed.  Since  the  analysis  becomes  increasingly  unstable  as  the 
mesh  is  refined,  it  is  necessary  to  maintain  as  coarse  a  mesh  as  possible.  For 
a  spherical  nose,  a  5x4  y-y  mesh  yielded  very  good: -results.  In  general,  even 
for  the  more  complex  nose  geometries,  a  10  x  12  T-y  mesh  should  be  a 
rough  upper  limit  of  mesh  density.  The  program  cannot  accommodate  a 
mesh  more  refined  than  15  x  15. 

Columns  3-10  indicate  the  y  mesh  interval  adjacent  to  the  bow  shock 
(  y=  1  ),  and  cols.  43  -  50  refer  to  the  y  mesh  increment  adjacent  to  the  axis 
of  symmetry.  Refer  to  Fig.  II  and  the  section  on  sample  inputs  and  outputs 
for  further  details. 

Thus  the  number  of  input  cards  necessary  for  a  given  run  varies  from 
roughly  8  to  an  upper  limit  of  22;  in  most  cases  the  number  will  nctexceed  12. 


81 


The  option  of  using  varying  mesh  intervals  was  included  so  that  a 
relatively  dense  mesh  may  be  prescribed  in  those  portions  of  the  subsonic  region 
where  the  velocity  gradients  are  high,  without  impairing  the  stability  (and  accuracy) 
of  the  solution.  A  recommended  mesh  for  a  spherical  nose  is: 

4  intervals  @  ^y;  =  0.  15  and  4  @  Ay2  =  0.1 

4  intervals  @  =  0.  15  and  4  @  AT2  =0.1 

These  proportions  result  in  a  relatively  coarse  (8x8)  mesh  which 
is  sufficiently  detailed  in  those  areas  characterized  by  high  velocity  gradients. 

For  those  nose  geometries  which  exhibit  a  very  rapid  change  of 
curvature  in  the  neighborhood  of  the  sonic  point  (e„g.  the  Apollo  configuration), 
it  will  be  necessary  to  utilize  3  mesh  intervals  in  each  direction  to  obtain 
accurate  results. 

Note  that  the  last  "body  card"  (card  A)  should  be  blank.  Thus,  there 
may  be  up  to  12  body  profiles,  followed  by  a  blank  card,  followed  by  the  nose 
geometry,  etc.  The  value  of  xt  for  all  body  profiles  must  be  non-zero. 


82 


TRANSONIC- SUBSONIC  PROGRAM 


The  input  format  is  identical  in  all  respects  to  that  used  for  the 
Flow  Field  program.  Since  only  the  nose  region  is  considered,  it  is  not 
necessary  to  prescribe  the  body  profiles  for  the  after-body  (cards  4,  5,  .  .  .), 
since  the  program  will  not  utilize  this  information.  Thus  only  one  or  two  body 
cards  which  prescribe  the  body  profile  in  the  vicinity  of  the  assumed  sonic  point, 
will  be  necessary. 

Refer  to  the  description  of  input  format  for  the  Flow  Field  program  for 
all  necessary  details.  Since  this  object  deck  is  half  the  size  of  the  Flow  Field 
deck,  it  may  not  be  advantageous  to  utilize  the  option  to  read  data  cards,  on¬ 


line. 
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SUPERSONIC  PROGRAM 


Option  1 


Option  1  is  employed  when  it  is  necessary  to  punch  all  input  data 
on  cards.  This  situation  occurs  whenever  the  binary  tape,  generated  by  the 
Transonic  -  Subsonic  program,  is  not  available.  Before  generating  the 
supersonic  flow  field  mesh,  the  program  will  write  all  inputs  on  binary  tape 
mounted  on  unit  B3,  in  exactly  the  same  format  as  that  required  by  options 
2  and  3  of  this  program,  for  subsequent  use,  if  desired. 


CARD  COL. 

DATA 

FORMAT 

1  1 

0 

I 

2-72  - 

Hollerith  information  -  title 

2  1 

1 

I 

3  1-10 

M 

00 

F 

11-20 

p 

F 

00 

21-30 

T 

00 

F 

31-40 

(  S/R) 

F 

oo 

41-50 

€o 

F 

51-60 

p 

00 

E 

4, 5, ..A  1-10 

a 

F 

11-20 

b 

F 

21-30 

c 

F 
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CARD  COL. 

DATA 

FORMAT 

31-40 

d 

F 

41-50 

e 

F 

51-60 

x 

F 

t 

61-70 

a 

F 

A+l  1-72 

Blank 

A+2.A+3 . B  1-8 

X 

F 

9-16 

y 

F 

17-24 

p> 

F 

25-32 

6 

F 

33-40 

S/R 

F 

41-48 

h' 

F 

49-56 

P' 

F 

57-64 

y 

F 

65-72 

p 

F 

B+l  1-72 

Blank 

Card  1  is  an  identification  card  which  is  printed  as  output.  The 
option  is  prescribed  on  card  2.  Card  3  contains  the  free  stream  conditions  and 
the  angle  (in  radians)  of  inclination  with  respect  to  the  x  axis  of  the  bow  shock 
at  the  reference  line  (  ).  Cards  4  through  A  contain  the  geometrical  data 

of  sequential  sections  of  the  body  profile;  a  maximum  of  12  is  permitted. 
Properties  of  sequential  points  along  the  reference  line  from  detached  shock  to 
body  are  punched  in  cards  A+2  through  B.  While  a  maximum  of  21  points  will 
be  accepted,  no  more  than  a  dozen,  or  so,  is  necessary  in  most  cases. 
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Conical  Nose 


Option  1  also  includes  the  case  of  a  body  having  a  conical  nose. 

In  this  case,  the  program  computes  the  points  along  a  reference  line  having 
a  constant  value  of  y.  Data  cards  1  through  A  +  1  are  identical  to  those  previously 
described;  cards  A  +  2  through  B  +  1  do  not  apply.  Note  that  the  values  of  a  and 
b  on  card  4  must  be  zero  since  the  first  body  profile  is  conical. 

Since  this  reference  line  is  written  on  the  binary  tape  on  unit 
B  3,  subsequent  runs  may  be  executed  using  options  2  and/or  3. 

Note  that  ^  on  card  2  must  be  left  biank  since  it  i-s  not  known 

as  an  input. 
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Option  2 


Option  2  is  employed  when  the  binary  tape  written  by  the  Transonic  - 
Subsonic  program  is  available,  and  the  supersonic  portion  of  the  flow  field  has 
not  been  completely  generated.  This  situation  could  occur  if  the  Transonic  - 
Subsonic  program  were  run,  or  if  the  Flow  Field  program  were  manually 
terminated  at  the  completion  of  Link  1  due  to  scheduling  problems,  or  if  a 
machine  failure  occured  during  execution  of  Link  2,  or  any  similiar  cause. 

When  this  option  is  chosen,  the  program  will  behave  essentially  like  Link  2 
of  the  Flow  Field  program;  binary  tape  B3  will  be  read  and  the  supersonic 
flow  field  mesh  will  then  be  constructed. 


CARD  COL.  DATA 

1  1  0 

2-72  Hollerith  information  -  title 

2  1  2 


FORMAT 

I 


I 


Card  1  may  be  used  for  identification  since  it  is  printed  as  output. 
The  option  is  prescribed  on  card  2. 
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Option  3 


Option  3  is  employed  when  the  binary  tape  written  by  the  Transonic- 
Subsonic  program  is  available,  and  it  is  desired  to  alter  the  geometry  of  the  after¬ 
body  of  the  vehicle  under  consideration,  and  generate  the  corresponding  super¬ 
sonic  flow  field  for  an  entire  family  of  body  shapes  without  incurring  the  cost  of 
repeating  the  transonic- subsonic  calculations.  The  new  body  equations  will  be 
written  on  the  binary  tape  B-3  to  replace  the  old  ones. 

CARD  COLS,  DATA  FORMAT 

11  0  I 


2-72  Hollerith  information  -  title 


2  1 

3 

I 

3,4, ...A  1-10 

a 

F 

11-20 

b 

F 

21-30 

c 

F 

31-40 

d 

F 

41-50 

e 

F 

51-60 

X 

t 

F 

61-70 

a 

F 

A  +  1  1-72 

Blank 

Card  1  may  be  used  for  identificatication  since  it  is  printed  as  output 

The  option  is  pre  scribed  on 

card  2.  Cards 

3  through  A  contain  the  geometrical 

data  of  sequential  sections  of  the  body  profile;  a  maximum  of  12  is  permitted. 


OPERATING  INSTRUCTIONS 
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These  computer  programs  are  designed  for  use  with  the  standard 
IBM  FORTRAN  Monitor  Systems  which  are  commonly  used  with  the  IBM  709/7090/7094 
digital  computers.  Since  the  FLOW  FIELD  program  is  a  CHAIN  job  consisting  of 
two  LINKS,  any  variation  of  the  above  systems  which  does  not  accommodate 
CHAIN  jobs,  cannot  be  utilized.  The  TRANSONIC- SUBSONIC  and  SUPERSONIC 
programs  are  not  CHAIN  jobs  and  thus  are  not  subject  to  the  above  restriction. 

An  IOU  table  subroutine  is  included  as  part  of  the  object  decks; 
this  avoids  those  difficulties  associated  with  varying  logical  tape  assignments 
from  one  computer  installation  to  another.  The  following  tape  units  are  common 
to  all  programs: 

A  1  IBM  FORTRAN  System  Tape 

A  2  Program  INPUT  Tape 

A3  Program  OUTPUT  Tape  -  to  be  listed. 

A 4  Intermediate  Monitor  tape  -  used  to  stack  LINKS 

of  CHAIN  job  for  FLOW  FIELD  program 
B  1  Intermediate  Monitor  Tape 

B  2  Intermediate  Monitor  Tape 

B  3  Binary  Tape  utilized  by  all  these  programs 

B4  "Clock"  Tape  -  utilized  by  Monitor  System 
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FLOW  FIELD  and  TRANSONIC  -  SUBSONIC  PROGRAMS 
Input  Option  No.  1 

If  desired,  the  object  decks  may  be  stacked  with  those  of  other 
programs  and  loaded  onto  tape  via  the  IBM  1401  as  part  of  a  "Monitor  Run".  In 
this  case,  the  data  (input)  cards  are  placed  immediately  behind  the  DATA  control 
card  and  become  a  part  of  the  object  deck  for  a  given  run.  Of  course,  this 
program  alone  may  be  loaded  onto  tape  as  indicated  above,  or  via  the  on-line 
card  reader,  if  desired.  The  tape  is  mounted  on  unit  A 2  as  the  INPUT  tape. 

This  procedure  is  standard  at  most  computer  installations;  the  lone  disadvantage 
is  that  the  binary  program  cards  must  be  handled  each  time  a  run  is  to  be 
executed. 

The  following  sense  switch  setting  is  required  for  this  input  option: 

Sense  Switch  5:  UP 

Input  Option  No,  2 

Since  the  object  deck  of  the  FLOW  FIELD  program  is  rather  large 
and  thus  somewhat  cumbersome  to  handle,  the  user  may  elect  to  store  this  deck 
permanently  on  magnetic  tape.  This  tape  is  then  mounted  on  unit  A2  as  the 
INPUT  tape,  and  the  following  simple  procedure  executed: 

1.  Clear  memory  and  read  in  via  the  on-line  card  reader,  the 
FORTRAN  START  card. 

2.  After  the  Monitor  starts  functioning  (a  few  seconds),  FEED  out 


this  START  card. 
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3.  Place  the  data  cards  in  the  hopper  of  on-line  card  reader  and 
depress  the  reader  START  key  until  the  READY  light  goes  on. 

The  following  sense  switch  setting  is  required  for  this  input  option: 

Sense  Switch  5:  DOWN 

Output  Options: 

Since  writing  BCD  tape  is  relatively  costly,  the  program  is  designed 
to  restrict  the  amount  of  output  and  generate  only  the  final  results  of  the  calculations. 
In  some  instances,  the  user  may  elect  to  generate  additional,  intermediate,  infor¬ 
mation  so  that  he  may  follow  the  development  of  the  computations  in  greater  detail. 
These  options  pertain  only  to  the  TRANSONIC  -  SUBSONIC  portion  (LINK  1)  of  the 
FLOW  FIELD  program;  the  SUPERSONIC  portion  (LINK  2)  makes  no  reference  to 
sense  switch  settings. 

Output  Option  No.  1 

The  program  will  write  on  output  tape  A  3  only  the  final  results  of  the 
calculations.  These  include  the  calculations  of  the  final,  iterated,  shock  polynomial, 
the  final  reference  line,  and  the  final  results  of  the  properties  within  the  subsonic 
(elliptic)  region  of  the  flow  field.  In  addition,  the  program  will  generate  a  binary 
tape  on  unit  B  3  which  is  used  as  input  media  by  LINK  2  of  the  FLOW  FIELD 
program.  This  tape  may  also  be  saved  for  subsequent  use  by  the  SUPERSONIC 
program.  The  following  sense  switch  settings  are  required  for  this  output  option: 
Sense  Switch  1:  UP 


Sense  Switch  2:  UP 


91 


Output  Option  No.  2 

In  addition  to  the  output  described  for  option  no.  1  the  program  will 
write  on  output  tape  A3,  the  properties  behind  the  bow  shock  at  point*,  in:  the 
t  ransonic  region,  the  properties  of  points  along  the  reference  line,  and  the  properties 
within  the  elliptic  region,  for  each  sweep  through  the  transonic- subsonic  regions 
in  the  process  of  iterating  to  obtain  the  final  shock  polynomial.  In  addition,  the 
coefficients  of  the  ba$ic,  perturbed,  and  final  shock  polynomials  are  written,  as 
well  as  statements  identifying  envelopes  in  the  transonic  region  (if  any).  The 
following  sense  switch  settings  are  required  for  this  output  option: 

Sense  Switch  1:  UP 

Sense  Switch  2:  DOWN 

Output  Option  No.  3 

In  addition  to  the  output  described  for  both  options  1  and  2,  the  program 
will  write  on  output  tape  A3,  the  properties  at  points  within  the  transonic  region 
for  each  sweep  through  the  transonic- subsonic  region.  This  writing  consumes  a 
substantial  amount  of  machine  time;  thus  this  option  should  be  utilized  very  sparingly. 
The  following  sense  ^witch  aettings  are  required  for  this  output  option: 

Sense  Switch  1:  DOWN 

Senae  Switch  2:  DOWN 

Any  combination  of  sense  switch  settings  involving  switches  1,2,  or  5 
may  be  selected  according  to  the  needs  of  a  particular  run.  Note  that  sense 
switches  3.4,  and  6  are  not  interrogated  by  the  program. 
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SUPERSONIC  PROGRAM 


As  indicated  above,  this  program  does  not  interrogate  sense  switch 
settings,  and  does  not  incorporate  any  input  and/or  output  options  as  is  present 
for  the  other  two  programs.  This  program  is  designed  to  function  as  a  monitor 
job;  the  data  cards  are  placed  behind  the  DATA  control  card  and  become  a  part  of 
the  object  deck  for  a  given  run.  This  deck  is  then  loaded  onto  magnetic  tape  via 
the  IBM  1401  computer,  or  via  the  on-line  card  reader  -  this  may  also  be 
accomplished  as  part  of  a  larger  monitor  run. 

There  are  three  variations  in  the  use  of  this  program;  these  are 
described  in  detail  in  the  section  on  input  formats  and  are  classified  as  execution 
options.  For  options  2  and  3,  the  program  requires  a  binary  input  tape  mounted 
on  unit  B  3;  this  tape  must  have  been  previously  generated  by  either  the  FLOW  FIELD 
program  or  the  TRANSONIC-SUBSONIC  program  or  generated  by  a.previous  run 
of  the  SUPERSONIC  program  -  see  below.  Option  number  1  requires  that  a  blank 
tape  be  mounted  on  unit  B  3  for  binary  output  purposes.  The  user  must  indicate 
whether  he  wishes  to  retain  the  use  of  the  tape  on  unit  B  3  for  future  runs  with  the 
SUPERSONIC  program. 


/7<r  I  'SCti£MAT/C  OF  now  F~/£LD 
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NOMENCLATURE 


This  section  will  define  all  symbols  and  their  units  which  have  not 
been -adequately  defined  within  the  context  of  the  description  of  the  analysis. 

Since  some  symbols  have  been  utilized  in  both  the  Transonic-Subsonic 
analysis  as  well  as  the  Supersonic  analysis,  and  may  have  different  connotations, 
this  table  of  nomenclature  is  subdivided  as  indicated  to  avoid  confusion., 


Transonic- Subsonic  Analysis 


P 

y 

6c 

u 

n 


T 

4» 


c 

U 

V 


Density,  slugs/cu.ft. 

Adiabatic  exponent 

Detachment  distance  at  axis 

Velocity  component  in  x  direction,  ft.  /sec. 

Dimensionless  pressure 

Radius  of  curvature  of  body,  at  axis 

Radius  of  curvature  of  detached  shock,  at  axis 

Velocity  component  in  y  direction,  ft.  /sec. 

(Basic  Shock  Folynominal  Analysis) 

TiP'y 

Dimensionless  stagnation  enthalpy 
Dimensionless  velocity  component  in  x  direction 
Dimensionless  velcoity  component  in  y  direction 


98 


0  Velocity  direction  measured  from  x  axis,  radians 

jU  Mach  angle,  radians 

f  (  T  ,  y  )  =  F  (output  listing)  Dimensionless  function  describing  entropy 

distribution 

P  Dimensionless  pressure 


R  =  P  =  6 


Dimensionless  density,  6 


y-  1 
y  +  1 


Subscripts: 

w 

o 

b 

Superscript: 

* 


Property  behind  detached  shock 
Property  at:  axis  of  symmetry 
Property  at  body 

Property  at  assumed  sonic  point 


Supersonic  Analysis 

W  Resultant  Velocity,  ft. /sec. 

M  Mach  Number 

y 

a 

Kj/2 


Adiabatic  exponent 
Speed  of  sound,  ft. /sec. 
Stagnation  enthalpy,  ft.  /  sec. 
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T 

P 

P' 

P 

I 

P 

h 

h' 


u 

v 

€ 


m 

S, 

'R 

Aij ,  Bmn 
0 
M 


Subscript: 


Spe  cific  heat  at  constant  pressure 

Reference  temperature  for  vibrational  excitation  of  air,  °R 

Stagnation  temperature  including  effects  of  vibrational  excita¬ 
tion  ,  °R 

Absolute  temperature,  °R 
Density,  slugs/cu.ft. 

Dimensionless  density 
Pressure,  Ibs/sq.ft. 

Dimensionless  pressure 
Enthalpy,  ft  /sec 
Dimensionless  enthalpy 

Velocity  component  in  x  direction,  ft.  /  sec. 

Velocity  component  in  y  direction,  ft.  /sec. 

Angle  of  inclination  of  shock  with  respect  to  x  axis,  radians 
Mass  flow  slugs/sec. 

Entropy 

Coefficients  of  numerical  fit  of  the  Mollier  diagram  for  air 
Velocity  direction  measured  from  x  axis,  radians 
Mach  angle,  radians 


oo 


Free  stream  condition 
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APPENDIX  I 


CHECK  OF  MASS  FLOW  ALONG  REFERENCE  LINE 


For  the  purpose  of  establishing  an  accurate  relationship  between  mass  flow 

and  entropy,  it  is  necessary  to  integrate  mass  flow  through  the  shock  layer  along  the 

reference  line,  from  body  (point  p)  to  detached  shock  (  point  1)  .  This  integrated 

value,  mj  ,  is  compared  with  the  mass  flow  crossing  the  detached  shock  between 

the  axis  and  point  1:  mj  =  V2  yZj  .  The  error  has  been  found  to  be  very 

small,  but  for  purposes  of  consistency  it  is  "carried  along"  in  all  subsequent 

calculations  of  mass  flow  across  the  detached  shock.  Thus  we  define  this  error, 
m,  # 

C  =  m  -  1:  then  the  mass  flow  crossing  the  detached  shock  at  any  point,  N, 

is  calculated  as  m  =  V,  p  W  (  y2  +  €  y2  ) 

N  00  co  N  . 


The  numerical  integration  of  mass  flow  along  the  reference  line  is 
calculated  as  follows: 


m 


1 

E7* 

j=P-l 


j+1  j+1 

Pj  W.  ain  (  Qj  -  )yj  +  #J+1Wj+1  sin  (  ftJ+r  <p.  )  yj+1 


(  W*  +  (yi 


“W 


2  1 


where 


j+1  -1 

<p .  =  tan 


7. 


-  x 


3+1 


x  . 
J 


APPENDIX  II 


FLOW  DIAGRAMS 


Most  copies  of  this  report  contain  only  4  pages  of  flow  diagrams 
which  describe,  in  general,  the  logical  flow  of  these  programs.  These  are 
contained  on  pages  FI,  F21,  F28  and  F49. 

A  few  copies  include  very  comprehensive  and  detailed  flow  diagrams 
of  these  programs.  These  are  included  for  reference  by  those  expert  program¬ 
mers  who  have  the  coding  listings  available,  and  who  wish  to  study  them  in 
some  detail. 

Those  who  wish  to  modify  the  coding  should  note  that  the 
TRANSONIC-SUBSONIC  program  (LINK  1)  consumes  virtually  all  of  core 
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APPENDIX  III 


SAMPLE  INPUT  AND  OUTPUT  LISTINGS 


Some  typical  sample  listings  of  input  data  cards  and  of  program  output 
are  included  and  described  below: 

Page  S  1  contains  a  listing  of  input  data  cards  to  the  SUBSONIC 
(and  FLOW  FIELD)  Program  at  the  top  of  the  page.  The  listing  at  the  bottom 
of  the  page  is  of  the  input  cards  to  the  SUPERSONIC  program,  option  number  1. 

Pages  S2  through  S5  contain  the  entire  output  of  the  SUBSONIC  Program, 
when  sense  switches  1  and  2  are  both  UP.  Page  S2  consists  of  the  input  data, 
which  is  written  for  purposes  of  identification  and  checking.  Page  S3  consists 
of  the.  final  reference  line  which  is  scaled  so  as  to  be  compatible  with  the 
prescribed  body  profile.  This  scaling  precludes  the  possibility  of  mesh  crossings 
in  region  A  of  the  SUPERSONIC  Program,  due  to  any  small  inconsistancy, 
between  the  body  point  on  the  reference  line  generated  by  the  Transonic  analysis, 
and  the  prescribed  body  profile.  Pages  S4  and  S5  exhibit  the  properties  within 
the  Subsonic  region  of  the  flow  field. 

Pages  S6  through  S20  contain  a  few  sheets  of  the  output  listing  of  the 
SUPERSONIC  program.  Page  S6,  which  is  the  first  page  of  output,  lists  the 
input  data  cards  for  identification  and  checking.  Pages  S7,  S  8,  and  S9  follow, 
and  are  the  first  points  generated  by  the  program  in  region  A.  The  final  first 
family  characteristic  line  originates  at  the  body,  which  is  the  last  input  point 
on  the  reference  line  (  x  =  .  33118,  y  =  .6804  on  page  S6),  continues  with  the 
first  IN  (interior)  point  on  page  S  10,  and  terminates  at  the  DS  (detached  shock) 
point  on  page  S  11.  Region  B  follows  with  the  BD  (body)  point  on  page  S  11,  and 


S 12. 


Further  on  in  region  B  an  expansion  corner  is  detected  as  illustrated 
on  page  S  13.  Still  further  on  in  region  B,  a  reentrant  corner  is  detected  and 
the  beginning  of  the  resulting  secondary  shock  and  surrounding  flow  field  is 
shown  on  pages  S  14  through  S  17.  A  typical  second  family  characteristic  line 
in  region  C,  originating  at  the  detached  shock,  crossing  a  secondary  shock,  and 
terminating  at  the  body,  is  shown  on  pages  S  18,  S  19,  and  S20. 
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